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An  improved  understanding  of  the  dynamic  strength  of  geologic 
solids  surrounding  an  underground  explosion  is  needed  for 
improving  calculations  of  the  teleseismic  amplitudes  from  these 
explosions.  Studies  of  one-dimensional  compression  and  shear  waves  in 
impacted  solids  have  the  potential  for  providing  information  that  can 
be  used  to  assess  and  extend  existing  material  models.  Of  particular 
interest  to  teleseismic  calculations  are:  direct  determination  of 
shear  strength,  determination  of  dynamic  frictional  characteristics 
along  material  surfaces,  and  determination  of  whether  large  amplitude 
shear  waves  can  propagate  in  the  materials  of  interest. 

We  chose  to  study  polycrystalline  salt,  a  material  of  interest 
for  teleseismic  calculations.  We  also  performed  tv;o  experiments  on 
Westerly  granite  to  get  an  indication  of  how  compression  and  shear 
waves  propagate  in  a  hard  rock. 

Compression  and  shear  wave  experiments  were  conducted  on  pressed 
polycrystalline  salt  obtained  from  the  Lawrence  Livermore  Laboratory 
(LLL).  The  uniaxial  strain  experiments  provided  results  (Figure  3) 
similar  to  previous  LLL  work.  Combined  compression  and  shear  wave 
experiments  were  then  conducted  to  obtain  longitudinal  and  shear 
particle  velocity  profiles.  The  longitudinal  profiles,  under  combined 
compression  and  shear  loading,  gave  profiles  (Figure  5)  similar  to 
uniaxial  strain  profiles,  indicating  negligible  effect  of  shear  loading 
on  the  longitudinal  stress-volume  (o^“^0  relation.  Despite  the  small 
porosity  (1%  to  1.5%),  pore-collapse  was  evident  in  the  cf^-V  relations 
(Figure  16) . 


The  observed  shear  wave  profiles  (Figure  6)  are  of  loiJ  amplitude 
and  shov/  strong  attenuation  with  propagation  distance.  Beyond  about 
7  mm,  the  shear  wave  amplitude  is  comparable  to  experimental  uncertain- 
ities.  The  start  of  the  shear  wave  signal  and  some  of  the  structure  in 
the  profiles  have  not  been  completely  explained  and  are  most  likely  an 
experimental  artifact  related  to  the  preceding  longitudinal  wave.  When 
these  experiments  were  repeated,  similar  results  were  obtained.  The 
main  conclusion  from  these  experiments  is  that  significant  shear 
amplitudes  cannot  be  propagated  into  salt. 

A  constitutive  model  was  constructed  for  salt  using  a  porous 
solid  description  for  the  pressure-volume  relation  and  an  isotropic 
strain-hardening,  viscoplastic  description  for  deviator  stresses.  The 
static  model  parameters  are  completely  determined  from  the  static  data 
(Figures  12  and  13).  The  strength  model  parameters,  determined  from 
uniaxial  stress  data,  give  good  predictions  for  uniaxial  strain  loading 
(Figure  14).  Dynamic  calculations  were  made  using  a  single  relaxation 
time  in  the  viscoplastic  model.  The  relaxation  time  was  determined  by 
matching  wave  profiles  from  one  uniaxial  strain  experiment.  Subsequent 
calculations  with  the  same  set  of  parameters  gave  good  agreement  with 
longitudinal  profiles  (Figure  17  and  18)  at  other  stress  and  under 
combined  compression  and  shear  loading.  The  calculated  shear  wave 
profiles  showed  rapid  attenuation  of  shear  wave  amplitude  and  strong 
localization  of  shear  strain  (Figure  18).  These  results  are  similar 
to  the  experimental  observations.  However,  the  structure  and  the 
arrival  time  do  not  match  the  experimentally  observed  profiles.  The 
reasons  for  this  mismatch  are  probably  experimental  uncertainties  as 
indicated  earlier. 

A  simple  analysis,  using  the  governing  equations,  reveals  that 
the  observed  localization  is  a  consequence  of  the  kinematic  conditions 
and  material  behavior.  The  localization  of  shear  strain  does  not 
permit  large  shear  wave  amplitudes  to  be  coupled  into  the  specimen 
interior. 


The  main  conclusions  from  the  salt  work  are:  the  static  measurements 
of  strength  are  adequate  for  use  in  dynamic  calculations;  the  small 
amiiiitude  and  the  rapid  attenuation  of  the  shear  wave,  though  interesting, 
limit  the  practical  usefulness  of  combined  compression  and  shear  wave 
studies. 

Of  the  two  experiments  on  Westerly  granite,  one  gave  good  results. 

The  shear  wave  profiles  indicate  a  large  shear  stress  amplitude  and  an 
arrival  time  that  agrees  with  extrapolation  of  ultrasonic  data.  The 
granite  results  suggest  that  compression  and  shear  traves  can  be  used 
for  studying  hard  rocks  to  determine  dynamic  strength  and  friction 
behavior. 

Based  on  the  work  done  to  date,  it  is  recommended  that  future 
shear  wave  studies  be  performed  on  hard  rocks.  These  studies  can 
address  the  following  objectives:  determination  of  the  dynamic  mean 
stress-volume  relation  and  determination  of  dynamic  friction  properties 
across  interfaces.  Measurements  of  shear  particle  velocity  amplitudes 
across  interfaces,  at  different  stress  levels,  should  be  helpful  in 
determining  the  role  of  friction  under  inelastic  deformation. 
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1.  INTRODUCTION 


BACKGROUND 


The  teleseismic  amplitudes  from  a  given  underground  explosion  depend 

1  2 

on  the  near  source  material  properties.  ’  Thus,  numerical  calculations 
of  the  teleseismic  radiation  field  require  an  accurate  description  of 
the  material  behavior  at  high  stresses  and  strain  rates.  In  particular, 
the  dynamic  strength  of  the  material  is  important  in  defining  the  extent 
of  the  inelastic  region  and  the  attenuation  through  this  region.  An 
improved  understanding  of  the  dynamic  strength  of  geologic  solids  is 
the  objective  of  the  present  work. 


Constitutive  models  used  to  describe  the  behavior  of  near  source 
material  are  based  largely  on  laboratory  data  from  ultrasonic,  quasi¬ 
static,  and  shock  wave  experiments.  Ultrasonic  data  provide  characteriza¬ 
tion  in  the  elastic  range.  Quasi-static  tests  (generally  triaxial 
compression)  characterize  the  material  in  the  inelastic  range.  These 
tests,  conducted  for  varying  loading  conditions,  serve  as  the  main 
input  to  the  strength  models.  Shock  wave  data  are  used  to  assess  the 

importance  of  high  loadjng  rates  and  provide  data  at  stresses  inaccessible 

3 

by  quasi-static  tests.  Shock  data  have  two  main  restrictions:  the 

stress  tensor  is  Incompletely  determined,  and  the  loading  is  restricted 
to  a  particular  path  in  strain  space.  Because  of  these  restrictions, 
shock  wave  data  are  used  indirectly;  that  is,  material  models  are 
iteratively  fitted  to  the  wave  profiles. 


The  loading  path  followed  by  material  near  an  underground  explosion 

is  not  the  same  as  the  loading  path  followed  in  the  laboratory  tests. 

Hence  data  obtained  under  more  general  loading  paths,  particularly  under 

dynamic  loading  conditions,  are  desirable  for  assessing  and  extending 

existing  material  models.  One  method  for  obtaining  material  response 

under  more  general  loading  conditions  is  the  study  of  one-dimensional 

4 

compression  and  shear  waves  in  impacted  solids.  These  studies  have  the 


11 
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potential  for  providing  several  new  pieces  of  information  that  are 
important  to  calculations  of  underground  explosions:  direct  determination 
of  shear  strength  to  improve  descriptions  of  dynamic  strength  and  inelastic 
deformation  in  current  material  models,  determination  of  dynamic  frictional 
characteristics  along  material  surfaces,  and  determination  of  whether 
large  amplitude  shear  waves  can  propagate  in  the  materials  of  interest. 

The  last  two  measurements  are  particularly  relevant  to  ascertaining  the 
effects  of  faults  and  joints  on  wave  propagation. 

SCOPE  OF  THE  WORK 

Propagation  of  large  amplitude  compression  and  shear  waves  in  solids 
is  a  new  subject  with  no  existing  data  on  materials  of  Interest.  Hence 
the  present  work  is  of  an  exploratory  nature.  The  emphasis  is  on  data 
acquisition  and  interpreting  the  data  in  terms  of  a  constitutive  model. 

We  chose  to  study  polycrystalline  salt,  a  material  of  interest  to  the 
teleseismic  work.  The  specific  objectives  of  the  work  were  to 

•  Determine  how  compression  and  shear  waves  propagate  in 
impact-loaded  salt. 

•  Develop  a  constitutive  model  for  salt  based  on  existing 
data  and  check  the  model  calculations  against  the  experi¬ 
mentally  observed  compression  and  shear  wave  profiles. 

Another  material  of  interest  to  the  teleseismic  work  is  granite.  To  get 

an  indication  of  how  compression  and  shear  waves  propagated  in  a  hard 

rock,  we  attempted  two  preliminary  experiments  on  Westerly  granite. 


2 .  EXPERIMENTS 


Six  plate  impact  experiments  were  performed  on  pressed  polycrystalline 

salt  (NaCl)  and  two  experiments  on  Westerly  granite.  The  experimental 

method  and  results  of  these  experiments  are  described  here.  All  the 

experiments  were  performed  using  the  IMPS  method,  described  in  a  previous 
4 

report  and  summarized  in  Appendix  A. 

SALT  EXPERIMENTS 

We  chose  to  work  on  the  finely  grained  polycrystalline  pressed 

5  6 

salt  used  in  previous  static  and  shock  studies  performed  at  Lawrence 
Livermore  Laboratory  (LLL)  because  the  study  of  compression  and  shear 
waves,  reported  here,  is  a  new  area  of  study.  By  studying  a  material 
that  had  been  studied  in  the  past  by  more  conventional  methods,  x^e  could 
draxi;  on  the  past  work  to  determine  the  implications  of  the  new  data. 

From  the  data  acquired,  we  could  then  define  future  experiments  on  the 
natural  salt  of  interest  to  the  teleseismic  xjork. 

The  salt  used  in  our  x^ork  was  taken  from  a  pressed  cylinder  of 
12-inch  diameter  by  12  inches  high  (30.5  by  30.5  cm)  supplied  to  us  by 
Dr.  D.  B.  Larson  (LLL) .  Density  and  sound  speed  measurements  on  two 
pieces  cored  from  the  cylinder  gave  the  folloi<ring  results; 

3 

Initial  density  (p^)  =  2.135  g/cm 

Longitudinal  velocity  (Cj^)  =  4.48  mm/ps 

Shear  velocity  (C,)  =  2.64  mm/ps 

Within  the  experimental  uncertainty  (2%),  these  values  are  in  good 
agreement  x.?ith  the  values  reported  in  the  LLL  static  material  property 
xrork.  ^ 

Experimental  Method 

The  specimens  used  for  the  impact  experiments  were  taken  from 
cylindrical  pieces  (6.35  cm  diameter  by  5  cm  high)  that  were  cored  from 
the  large  cylindrical  block.  Because  of  the  hygroscopic  natxire  of  salt. 
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all  pieces  were  stored  in  a  desiccator  (or  vacuum)  when  not  in  use. 

The  cylindrical  pieces  x^ere  cut  into  thinner  discs  (betx^eeh  2  and  4  mm- 
thick)  and  were  ground  and  lapped  to  provide  flat  arid  parallel  surfacsri. 
Tolerances  were  typically  +0.05  mm. 

Preparation  oi  the  thin  discs  (less  than  3i5  mm)  was  mote  difficult 
tlian  anticipated  because  the  specimens  tended  to  develop  partial  cracks 
through  the  specimen  thickness.  Careful  gririding  and  lapping  procedures 
were  later  developed  to  obtain  the  thin  crack-free  specimens.  Thin 
specimens  are  necessary  for  performing  multiple  gage  experiments. 

The  specimen  assembly  shown  in  Figure  1  was  used  in  the  experiments. 
Several  thin  discs  (generally,  three  discs)  were  stacked  together  and 
backed  by  a  thicker  disc  (5  ram  or  greater) .  The  particle  velocity  gages 
consisted  of  thin  metallic  foils  (U-shaped)  with  ari  active  elemerit  at 
each  specimen  interface  and  gage  leads  coming  out  through  the  back.  A 
soft  aluminum  ribbon,  0.001  inch  thick,  was  bonded  into  an  accurately 
machined  groove  (0.040  inch  wide  and  0.001  Inch  deep)  pri  the  specimen 
surface  and  pulled  through  the  vertical  holes  (0.04  inch  in  diameter); 

The  holes  were  then  filled  with  a  mixture  of  salt  arid  epoxy  to  minimize 
perturbation  in  the  propagating  waves.  After  the  gage  was  emplaced,  the 
disc  surfaces  were  Lapped  to  ensure  that  the  gage  surface  was  flush  with 
the  specimen  surface.  The  complete  specimen  assembly  was  fabricated  fay 
sequentially  bonding  the  specimen  discs  using  an  epOxy  resiri;  The  weak 

* 

cohesive  strength  of  the  pressed  salt  permits  the  use  of  ari  epoxy  bond; 

From  the  brief  description  presented  above,  it  can  be  seen  that 
specimen  assembly  for  these  experiments  was  difficult  and  time  corisumirig; 
Hoii’ever,  this  procedure  xv’as  necessary  to  ensure  that  the  gage  was 
surrounded  by  the  salt  and  to  avoid  gage  lead  artifacts. 

After  bonding,  the  assembled  discs  were  potted  in  a  target  ring 
and  impacted  by  either  a  PbPlA  or  a  salt  flyer  plate.  A  schematic  view 
of  the  impact  experiment  is  shown  in  Figure  2.  Experimental  measurements 
consisted  of  monitoring  particle  velocity  time  profiles  at  various  depths 
in  the  specimen  plate.  In  additiorii  projectile  velocity,  impact  tiltj 
and  magnetic  field  strengths  were  monitored  in  each  e^eriment;  The 

To  avoid  slippage  at  the  interface,  the  interface  adhesion  shduld  Bi- 
higher  than  material  cohesion. 
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(a)  SIDE  VIEW  (b)  TOP  VIEW 

MA-5746-48 

FIGURE  1  SALT  SPECIMEN  ASSEMBLY  SHOWING  THE  GAGE  ELEMENTS  AND  HOLES 
FOR  LEADS 

The  three  thin  and  thick  backing  plates  are  shown  in  the  side  view. 


resu.1  of  the  six  impact  experiments  conducted  on  salt  are  discussed 
next  These  results  are  analyzed  and  compared  with  material  model 
predictions  in  Section  4. 


Table  1  lists  details  of  the  six  experiments.  Two  types  of  experi¬ 
ments  were  performed:  uniaxial  strain  (US)  experiments  and  experiments 
under  combined  compression  and  shear  (C  +  S) .  For  the  combined  compression 
and  shear  experiments,  the  impact  surfaces  were  inclined  at  75°  to  the 
projectile  axis  (see  Figure  2) .  In  contrast,  this  angle  was  90°  for 
uniaxial  strain. 

Results  of  Uniaxial  Strain  Experiments 

Experiment  1  (78-2-25)  was  a  uniaxial  strain  experiment  with  a  PMMA 
flyer  plate  and  an  impact  velocity  of  0.407  mm/ps.  The  magnetic  field 
was  aligned  in  the  usual  manner  to  measure  the  longitudinal  particle 
velocity,  u^(t).  This  experiment  was  similar  to  an  LLL  experiment. 

Figure  3  shows  the  voltage-time  profiles  from  the  five  gages.  Because 
the  voltage  and  time  scales  are  the  same  in  all  pictures,  direct  comparisons 
can  be  easily  made.  The  gage  locations  were  gage  1  (0  mm),  gage  2 
(3.79  mm),  gage  3  (7.52  mm),  gages  4  and  5  (11.22  mm).  The  leads  from 
gages  1  through  4  came  out  of  the  side  and  the  leads  from  gage  5  came 
out  through  the  back.  The  side  leads  cause  a  small  increase  in  signal 
before  the  arrival  of  the  release  wave,  as  seen  in  the  first  four  records. 
Gages  4  and  5  show  identical  profiles  before  the  small  increase  in  gage  4. 

The  first  gage  shows  a  flat  top  signal  (except  for  a  minor  perturbation 
* 

at  the  start  )  expected  from  the  impact  of  PMMA  on  salt.  The  subsequent 
gages  show  a  two-wave  structure:  a  small  elastic  precursor  followed  by 
a  large  plastic  wave.  These  profiles  are  in  good  agreement  with  the 
LLL  work.^  The  following  measurements  were  obtained  from  this  experiment: 
Peak  particle  velocity  =  0.14  +  0.01  mra/ps 

Elastic  wave  velocity  =  4.33+0-1  mm/ps 

Release  wave  velocity  =  4.94+0.15  mm/ps 

* 

This  occurs  if  the  gage  is  slightly  above  the  surface  of  the  specimen. 

*1* 

Normally  this  error  is  2%  to  3%.  The  higher  value  quoted  reflects  some 
calibration  difficulties  in  this  particular  experiment. 
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The  first  two  values  are  in  good  agreement  with  the  LLL  irork.  The 

release  wave  velocity  (in  Lagrangian  coordinates)  is  significantly 

0 

different  from  the  LLL  value  of  5.73  mm/ps.  Because  of  the  procedure 

ic 

used  in  the  LLL  experiments,  we  believe  that  our  value  is  more  accurate. 

Our  value  also  agrees  better  with  ultrasonic  measurements  in  polycrystalline 
salt.^  Further  analysis  of  the  profiles  shown  in  Figure  3  are  presented 
in  Section  4. 

Experiment  2  (78-2-33)  was  also  a  uniaxial  strain  experiment,  but 
the  magnetic  field  was  parallel  to  the  direction  of  the  projectile  motion. 
The  purpose  of  this  experiment  was  to  determine  the  accuracy  of  our  shear 
wave  measurements  described  later.  In  a  uniaxial  strain  experiment,  only 
the  Uj^  component  of  particle  velocity  is  non-zero.  Hence,  if  the  magnetic 
field  is  aligned  parallel  to  u^^,  we  expect  a  zero  signal  under  ideal 
conditions.  However,  because  of  small  errors  in  tilt,  gage  emplacement, 
and  field  alignment,  we  can  expect  a  small  error  signal.  In  Appendix  A, 
we  have  estimated  the  magnitude  of  this  signal  to  be  2  percent  of  the 
longitudinal  signal. 

To  experimentally  determine  the  magnitude  of  this  error  signal, 
we  performed  experiment  2,  which  was  similar  to  our  C  +  S  shots  in  all 
respects  except  that  it  was  a  uniaxial  strain  experiment.  The  magnetic 
field  strength  used  was  the  highest  possible  value  so  that  it  would 
generate  as  large  a  signal  magnitude  as  possible.  Gages  were  located  at 
the  impact  surface  and  at  a  depth  of  8.02  mm.  The  gage  leads  in  this 
experiment,  as  in  the  C  +  S  shots,  came  out  through  the  back  of  the 
specimen.  Upon  exit,  these  leads  are  folded  to  make  electrical  connec¬ 
tion-.  Hence  the  useful  recording  time  of  the  gage  is  determined  by 
the  arrival  of  the  wave  at  the  rear  surface.  Data  beyond  this  time  are 
not  pertinent  to  the  results  and  have  not  been  shown  in  any  of  the 
photographs . 


In  the  LLL  work,  the  release  wave  comes  from  the  back  of  the  specimen 
and  increases  the  particle  velocity.  In  our  work,  the  release  wave 
comes  from  the  back  of  the  impactor  and  decreases  the  particle  velocity. 
The  decrease  in  particle  velocity  is  more  accurately  observed  than  the 
increase  because  lead  effects  can  also  give  rise  to  a  small  increase. 
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Figure  4  shows  the  particle  velocity  records  from  the  two  gages. 

Taking  into  account  the  baseline  shifts,  the  impact  surface  gage  shows 
a  +  13-mV  signal  and  the  interior  gage  shows  a  +  18-mV  signal.  These 
values  are  comparable  to  the  11-mV  signal  estimated  in  Appendix  A.  The 
observed  signal  corresponds  to  a  shear  particle  velocity  magnitude  that 
is  2%  to  3%  of  the  longitudinal  particle  velocity  value.  Hence  2%  to 
3%  of  the  longitudinal  particle  velocity  value  represents  the  accuracy 
limit  on  the  shear  velocity  aim)litudes  reported  here.  With  additional 
work,  the  field  alignment,  tilt,  and  gage  placement  can  be  improved  to 
reduce  this  error  to  less  than  a  percent. 

Results  of  Combined  Compression  and  Shear  Experiments 

Experiments  3  (79-2-27)  and  4  (78-2-28)  were  performed  with  the 
impacting  surfaces  inclined  at  an  angle  of  75“  to  the  projectile  axis. 

Both  experiments  had  salt  flyer  plates  and  an  impact  velocity  of  0.41  mm/ps. 
Each  experiment  had  4  gages  with  leads  coming  out  of  the  back.  For  the 
specimen  thicknesses  used,  a  useful  recording  time  of  3  ps  after  impa'-t 
is  obtained  in  these  experiments.  At  3  ps,  the  longitudinal  wave  reaches 
the  back  surface  and  the  leads  are  perturbed. 

In  experiment  3  (78-2-27),  the  magnetic  field  was  aligned  to  measure 
only  the  longitudinal  component.  The  gage  records  at  0,  2.56,  5.17,  and 
7.62  mm  are  shown  in  Figure  5.  A  step  input  at  the  impact  surface, 
resulting  in  a  two-wavs  structure,  can  be  seen  in  the  specimen  interior. 

The  two-wave  structure  is  less  distinct  than  that  observed  in  experiment  1 
(Figure  3),  and  the  plastic  wave  at  this  peak  stress  is  much  steeper  than 
the  plastic  wave  in  experiment  1.  Gage  3  in  Figure  3  and  gage  4  in 
Figure  5  are  nearly  at  the  same  location,  but  differences  in  the  wave 
profiles  can  be  observed.  The  less  distinct  separation  between  the  two 
waves  and  the  steeper  plastic  wave  suggest  that  the  relation  steepens 

significantly  in  the  stress  range  between  these  shots.  The  trace  beyond 
the  center  of  the  photograph  is  not  pertinent  because  of  lead  effects 
discussed  earlier. 


In  experiment  4  (78-2-28),  the  magnetic  field  was  aligned  to  measure 
only  the  shear  component.  The  gage  records  at  0,  2.12,  4.72,  and  7.31  ram 
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FIGURE  4  VOLTAGE-TIME  PROFILES  AT  VARIOUS 
GAGE  LOCATIONS  (EXPERIMENT  2. 
78-2-33,  UNAXIAL  STRAIN) 
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FIGURE  5 


VOLTAGE-TIME  PROFILES  AT  VARIOUS 
GAGE  LOCATIONS  (EXPERIMENT  3.  78-2-27 
COMBINED  COMPRESSION  AND  SHEAR) 


are  shown  in  Figure  6.  These  gage  locations  are  similar  to  the  gage 
locations  in  experiment  3.  Data  beyond  the  useful  recording  time  of 
the  gage  have  not  been  shown. 

The  first  gage,  located  at  the  ■>mpact  surface,  shows  the  particle 
velocity  amplitude  theoretically  expected  from  salt/salt  impact  and 
assuming  no  slip.  The  rise  time  is  slower  than  in  the  longitudinal 
measurements.  This  result  has  also  been  observed  at  the  impact  surface 
in  other  materials  and  suggests  that  the  same  gage  takes  longer  to  reach 
peak  amplitude  in  shear  than  in  compression.  After  the  peak  value  is 
reached,  there  is  some  oscillation  about  the  peak  amplitude.  Although 
we  do  not  have  a  explanation  for  this  phenomenon,  it  is  similar  to  the 
'stick'  and  'slip'  phenomenon  observed  in  frictional  studies.^ 

The  interior  gages  show  considerably  smaller  amplitudes  and  a 
complex  structure.  Because  of  the  small  amplitude  of  these  profiles, 
we  first  compare  the  amplitude  of  the  shear  signals  with  experimental 
errors  determined  in  experiment  2  (78-2-33).  To  eliminate  differences 
resulting  from  differences  in  magnetic  field  strengths  and  gage  lengths, 
we  make  these  comparisons  by  normalizing  the  shear  wave  amplitudes  with 
respect  to  the  longitudinal  particle  velocity.  The  first  peak  in 
particle  velocities,  shown  in  Figure  6,  represent  27%,  10%,  7.5%,  and 
5%,  respectively,  of  the  longitudinal  particle  velocity  imparted  at  the 
impact  surface.  Although  the  amplitude  measured  in  gage  4  is  not  much 
larger  than  the  2%  to  3%  error  observed  in  experiment  2  (78-2-33),  the 
remaining  amplitudes  are  outside  the  error  limit. 

To  compare  the  Lime  of  arrival  of  the  shear  signal  with  the  longi¬ 
tudinal  wave  and  the  ultrasonic  x^ave  velocities,  we  have  drawn  a  h-t 
plot  in  Figure  7.  The  gage  locations  in  the  two  experiments  (3  and  4) 
were  not  identical,  but  are  close  enougli  to  interpolate  for  accurate 
comparisons.  In  Figure  7,  we  show  lines  corresponding  to  the  arrival 
of  the  longitudinal  wave,  onset  of  the  longitudinal  plastic  wave,  and 

"k 

90%  of  the  longitudinal  peak  for  experiment  3.  For  experiment  4,  we 
show  tlie  onset  of  the  shear  signal  and  the  first  peak  following  this 
signal.  These  I'esults  give  a  longitudinal  elastic  wave  velocity  of 

k 

Because  of  rounding,  determination  of  the  arrival  time  of  the  peak 
longitudinal  wave  is  difficult  and  the  90%  point  is  used. 


4.36  mm/us,  in  good  agreement  with  the  4.33  mm/jjs  measured  in  experiment 
1.  The  shear  wave  arrival  obtained  from  Figure  6  has  a  value  of 
3.55  +  0.15  mm/ps  and  lies  between  the  onset  of  the  plastic  wave  and 
the  90";  longitudinal  peak.  This  value  of  shear  wave  velocity,  even 
after  correction  for  compression,  is  much  higher  than  the  ultrasonic 
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shear  wave  velocities  measured  under  static  high  pressure  experiments. 

It  appears,  therefore,  that  the  onset  of  the  shear  signal  in  Figure  5 
is  an  experimental  artifact.  Because  of  the  small  magnitude  of  the 
shear  signal,  it  is  not  possible  to  get  an  accurate  time  of  arrival  for 
the  shear  wave.  Furthermore,  the  structure  in  gage  records  2  and  3  in 
Figure  6  makes  it  difficult  to  link  a  particular  break  in  the  slope 
with  shear  wave  arrival.  Hence  the  main  conclusion  we  can  draw  from 
experiment  4  is  that  the  shear  wave  amplitude  attenuates  very  rapidly 
within  1-2  mm  of  the  impact  surface,  and  by  about  7  mm,  the  signal 
amplitude  is  negligible. 

Experiments  3  and  4  were  repeated  to  confirm  the  experimental 
results.  The  results  from  these  experiments  (5  and  6)  v^ere  generally 
similar  except  for  the  shear  wave  profile  seen  at  the  impact  surface. 

Figure  8  shows  the  shear  xi?ave  records  from  experiment  6.  The  impact 
surface  gage  shows  an  initial  jump  to  the  expected  value  and  subsequently 
decay  to  a  much  lower  value.  Two  possible  causes  for  this  difference 
are:  (1)  the  slip  conditions  at  the  impact  surface  are  different  or 

(2)  the  gage  in  experiment  6  was  slightly  below  the  surface  and  the 
amplitude  decays  because  of  the  rapid  shear  attenuation  near  the  impact 
surface.  At  present,  we  cannot  distinguish  betv/een  these  two  possibilities, 
rhe  interior  gages  show  records  generally  similar  to  those  sho\m  in 
Figure  6.  Again,  no  clear  arrival  can  be  noted  for  the  shear  wave. 

The  results  from  these  experiments  can  be  summarized  as  follows: 
the  longitudinal  wave  profiles  under  combined  compression  and  shear  are 
similar  to  the  uniaxial  strain  profiles;  the  shear  wave  profiles  show 
very  large  attenuation  in  amplitude  near  the  impact  surface,  and  very 
little  shear  motion  is  coupled  into  the  specimen  interior.  The  accurate 
determination  of  shear  wave  arrival  is  not  possible  from  the  present 
data.  The  rapid  shear  wave  attenuation  is  discussed  further  in  Section  4 
and  compared  with  the  material  model  predictions. 
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FIGURE  B 


VOLTAGE-TIME  PROFILES  AT  VARIOUS 
GAGE  LOCATIONS  (EXPERIMENT  6,  78-2-29, 
COMBINED  COMPRESSION  AND  SHEAR) 


GRANITE  EXPERIMENTS 


The  purpose  of  the  granite  experiments  was  to  determine  the  feasi¬ 
bility  of  propagating  and  measuring  large  amplitude  compression  and 
shear  waves  in  a  hard  rock.  Two  experiments,  with  different  gage 
emplacement  techniques,  were  performed.  Only  one  of  these  was  successful 
and  is  described  below  (Proj .  Vel.  =  O.Al  mm/psec.  Angle  6  =  12°). 

Experimental  Method 

The  experimental  assembly  for  the  granite  shot  was  similar  to  that 

for  the  salt  shots.  Particle  velocity  gages  x^ere  inlaid  in  the  material 

with  gage  leads  coming  out  of  the  back  as  shoim  in  Figure  1.  Copper 

wires  (0.002  inch  in  diameter)  instead  of  aluminum  foils  were  used  as 

gage  elements.  These  small  wires  were  connected  to  0. 03- inch-diameter 

copper  wires,  x^hich  served  as  gage  leads.  The  space  around  the  gage 

* 

element  and  gage  leads  was  filled  Xi7ith  aluminum  Devcon  to  minimize 
perturbation  in  the  x^ave  front.  The  granite  discs  xvere  not  bonded 
together  xvith  epoxy  like  the  salt  discs,  but  were  mechanically  pressed 
together  and  potted  in  the  aluminum  ring.  It  was  felt  that  the 
adhesion  strength  caused  by  friction  betX(?een  the  granite  pieces  would 
be  stronger  than  the  cohesive  and/or  dhesive  strength  of  the  epoxy. 

Because  x^?e  xrere  able  to  perform  only  one  experiment,  x^e  xi7anted  to 
observe  signals  from  both  the  longitudinal  and  shear  wave.  This  x,'as 
accomplished  by  using  a  magnetic  field  configuration  that  gives  rise 
to  a  signal  from  each  xvave.  The  longitudinal  and  shear  signals  are 
distinguished  by  the  arrival  time  of  the  txro  waves.  This  approach  to 
distinguishing  the  tx^o  xvaves  is  applicable  only  when  distinct  and  large 
signals  are  expected  from  each  wave,  unlike  the  salt  results.  As 
discussed  next,  the  signals  from  the  granite  experiment  x<rere  indeed 
large  enough  to  permit  this  distinction. 

Experimental  Results 

The  particle  velocity  profiles  from  the  three  gages  are  shoxTO  in 

Figure  9.  As  expected,  these  records  display  substantial  piezoelectric 

noise.  The  first  gage  shox^s  a  square  wave  input,  the  second  gage 

starts  to  shox^  a  wave  separation,  and  the  third  gage  clearly  shows  the 

two-x,'ave  structure:  compression  and  shear  wave. 

_ 

A  commercially  available  mixture  of  aluminum  and  epoxy. 
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Figure  10  shows  the  voltage-time  signal  from  gages  1  and  3  along 
with  the  signal  expected  from  the  granite  for  a  no-slip,  elastic  response. 
The  plots  in  Figure  10(a)  show  that  the  measured  signal  is  in  good 
agreement  with  the  expected  signal  for  a  symmetric  no-slip  impact;  that 
is,  the  shear  velocity  transmitted  to  the  specimen  is  one-half  the 
projectile  velocity  along  the  shear  direction.  The  plots  in  Figure  10(b) 
show  that  the  first  wave  amplitude  agrees  with  the  expected  longitudinal 
signal.  However,  the  second  wave  amplitude  is  slightly  smaller  and 
considerably  more  dispersive  than  that  expected  from  an  elastic  response. 
The  amplitude  of  the  shear  tv’ave  is  a  measure  of  the  dynamic  strength,  and 
the  wave  profile  shape  indicates  a  convex  shear  stress-strain  relation. 

The  wave  velocities  used  in  generating  the  elastic  profiles  in 

Figure  10  are:  a  longitudinal  vsave  velocity  of  5.35  +  0.25  rnm/ps 

measured  in  these  experiments  and  a  shear  wave  velocity  of  3.67  rara/ys 

4 

obtained  from  extrapolation  of  ultrasonic  data. 

Because  of  the  quality  the  records,  further  analyses  were  not 
attempted.  The  main  conclusions  from  the  granite  shot  are  that  compression 
and  shear  waves  can  be  propagated  and  measured  in  hard  rocks.  By  extending 
the  stress  range  to  higher  values,  we  can  determine  the  shear  stresses 
that  cause  shear  failure  and  find  out  about  the  post-failure  behavior 
of  the  solid.  In  addition,  the  measurement  of  shear  wave  velocity  in 
the  shocked  state  can  allow  direct  determination  of  the  mean  stress-volume 
relation  in  the  shocked  state,  as  discussed  in  Section  4.  These  measure¬ 
ments  in  hard  rocks  may  also  allow  us  to  obtain  new  information  regarding 
dynamic  inelastic  deformation  (see  Section  4) .  Future  effort,  however, 
should  focus  on  a  hard  rock  without  the  piezoelectric  noise  to  allow 
more  accurate  investigation  of  material  effects. 


FIGURE  10  COMPARISON  OF  SIGNAL  DETERMINED  FROM  LINEAR 
ELASTIC  CONSIDERATIONS  AND  MEASURED  SIGNAL 


3.  CONSTITUTIVE  MODEL  FOR  SALT 


An  important  objective  of  the  present  work  was  to  develop  a 

mechanical  constitutive  model  for  salt.  Static,  ultrasonic,  acd  plate 

5  6 

impact  uniaxial  strain  data  obtained  at  LLL  ’  provided  the  basis  for 
the  form  of  the  constitutive  relations  and  the  specific  material 
parameters.  These  constitutive  relations  were  used  to  calculate  wave 
profiles  under  combined  compression  and  shear  loading  for  comparison 
with  experimental  data  obtained  in  the  present  work.  These  comparisons 
are  presented  in  Section  4. 

This  section  presents  the  main  features  of  the  constitutive 
relations  and  the  determination  of  the  material  parameters.  Detailed 
derivation  of  the  equations  and  incorporation  into  the  wave  propagation 
code  are  given  in  Appendices  B  and  C. 

FORMULATION  OF  THE  MODEL 

In  choosing  the  constitutive  description,  we  concentrated  on  using 
tlie  simplest  continuum  model  that  would  incorporate  the  experimentally 
observeu  features.  The  main  elements  of  the  model  were  determined 
from  the  static  data  with  later  modifications  for  incorporating  plate 
impact  data.  This  approach  was  used  because  the  static  data  are  more 
complete  than  the  dynamic  data. 

During  deformation,  salt  displays  many  features  typical  of  metal 

plasticity:  elastic-ductile  plastic  deformation,  strain-hardening, 

5  * 

and  lack  of  dilatant  behavior.  Uliile  porosity  effects  need  to  be 
included  in  the  constitutive  relations,  the  small  amount  of  porosity 
permits  considerable  simplification.  The  effect  of  shear  on  the  pressure- 


The  salt  specimens  obtained  from  LLL  have  1%  to  1.5Z  porosity. 

‘We  use  'pressure'  to  denote  the  negative,  ci  mean  stress.  For  purely 
hydrostatic  loading,  'pressure'  i^ill  be  replaced  by  'hydrostatic- 
pressure.  ' 
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volume  relation  is  small  and  can  be  approximated  by  using  the  pressure- 
volume  relation  obtained  from  static  uniaxial  strain  data.  Shear- 
enhanced  compaction  need  not  be  explicity  incorporated  in  the  pressure- 
volume  relation  and  in  the  subsequent  calculations.  The  thermal  effects 
of  pore-collapse  and  yielding  are  also  ignored  in  the  present  work. 

Because  of  the  above  considerations,  the  pressure  and  deviator 
parts  of  the  constitutive  relation  can  be  uncoupled  and  determined 
separately. 


a. . 


6. .  +  o'. 

ij  ij 


(1) 


where  o  /3  is  the  mean  stress  or  pressure  and  a',  are  the  stress 
mm  ij 

deviators.  Forms  of  the  equations  describing  the  stress-strain 
relations  for  each  of  these  two  parts  are  presented  in  the  next  two 
subsections. 

Of  the  various  materials  of  interest  to  the  DARPA  teleseismic 
program,  salt  comes  closest  to  a  metal  in  terms  of  its  deformation 
behavior.  This  similarity  simplifies  the  constitutive  description. 


PRESSURE-VOLUME  REIJVTION 

Although  the  porosity  of  the  salt  being  studied  is  small  (1%  to 
1.5Z),  both  static  and  dynamic  data  show  pore-compaction  effects. 

Hence  the  pressure-volume  relation  should  incorporate  pore-collapse 
and  consolidation  to  the  solid  material.  Of  the  various  models  used 
to  describe  porous  solids,  we  chose  the  empirical  relation  presented 
in  Seaman's  work  on  radiation  deposition  in  porous  metaJs.  By 
neglecting  thermal  effects,  we  have  further  simplified  this  model  for 
our  purposes.  A  brief  description  of  Seaman's  model  is  presented  below. 
A  listing  of  the  porous  P-V  model,  subroutine  POROUSl,  and  a  flow  chart 
demonstrating  how  the  equations  (which  are  presented  later)  are  used 
in  the  numerical  calculations  are  given  in  Appendix  B  along  with  the 
procedure  for  inserting  POROUSl  into  the  COPS  code. 
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Figure  11  shows  the  pressure-volume  relationship  typically  observed 
for  a  porous  solid.  As  will  be  shown  later,  a  similar  form  is  observed 
for  salt.  To  describe  the  pressure-volume  relation  shown  in  Figure  11, 
we  consider  initial  loading  in  the  porous  state  (path  ABC),  unloading 
and  reloading  in  the  porous  state  (BD) ,  and  unloading  and  loading  in 
the  solid  state  (CF) . 

Loading  in  the  porous  state  (curve  ABC)  is  characterized  by  an 

:5r 

initial  elastic  response,  pore-compaction,  and  consolidation  to  the 
solid  material.  In  our  model,  we  assume  that  the  mean  stress-volume 
relation  of  the  porous  solid  after  consolidation  is  identical  to  that 
of  the  nonporous  material.  This  assumption  is  only  an  approximation 
but  is  reasonable  within  our  use  of  the  model.'  In  the  model,  curve 
ABC  in  Figure  11  is  defined  as  a  series  of  parabolic  segments.  A 
three-segment  approximation  to  curve  ABC  is  shown  in  Figure  11.  Each 
segment  is  defined  by  an  equation  of  the  form 


P  =  P.  + 
1 


(V  -  V.) 


4DP 


(V.,.  -  V.)' 
1+1  1 


(V  -  V^)(V  -  V.^^)  (2) 


where,  V  is  the  specific  volume  and  the  subscripts  i  and  i+1  denote 
quantities  at  the  beginning  and  end  of  a  segment,  respectively.  The 
second  term  on  the  right  hand  side  of  Equation  (2)  represents  a  linear 
P(V)  variation  between  and  3nd  the  quantity  DP  is  the  maximum 

deviation  from  this  linear  curve.  In  segments  1,  2,  and  3  of  Figure  11, 
DP  is  positive,  zero,  and  negative,  respectively. 

Equation  (2)  was  derived  by  assuming  a  P(V)  variation  of  the  form 


P  =  a  +  b(V  -  V^)  +  c(V  -  V^)(V  -  V.^j^) 


(3) 


Elastic  implies  no  pore-compaction. 

^If  thermoraechanical  effects  are  included,  the  P-V-E  surface  of  the 
porous  and  solid  material  will  not  be  coincident  in  Che  P-V  plane. 


for  a  segment.  The  constants  a  and  b  were  determined  by  requiring 
that  P(V.)  =  P.  and  P(V..,)  =  P.,,;  the  constant  c  was  obtained  by 
solving  for  the  location  of  the  extremum  of  the  quadratic  term  of 
Equation  (3).  This  is  the  procedure  used  by  Seaman  to  model  the  initial 
loading  portion  of  a  porous  material  P-V  relationship  and  is  described 
in  detail  in  reference  9.  By  proper  choice  of  segments,  the  changes  in 
slope  at  segment  joints  are  minimized. 

Unloading  from  a  partially  consolidated  state  (path  BD)  is  assumed 
to  be  elastic,  that  is,  no  change  of  porosity.  The  unloading  path  is 
given  by 

AP  =  K  (4) 


where 


AP  ar.d  AV  =  pressure  and  specific  volume  increments 

V  =  average  of  the  specific  volumes  at  the  beginning 

and  end  of  a  cycle 

K  =  bulk  modulus  between  the  initial  and  consolidated 

states  defined  as 


K  -  K 

K  =  K-  +  — -  (p  -  p,) 

I  P_-Pt  avg  "'I-' 


where 

and  K  =  bulk  moduli  in  the  initial  and  consolidated  states 
I  c 

p_  and  o  =  initial  and  consolidated  densities 
1  ■  c 

p  =  average  of  the  densities  at  the  beginning  and  end 

of  a  cycle. 

Reloading  occurs  elastically  along  the  same  path  until  the  original 
loading  curve  is  reached.  Curve  BD  in  Figure  11  represents  a  typical 
unloading/reloading  path  in  a  partially  consolidated  state. 

If  the  pressure  is  high  enough,  all  the  porosity  is  expected  to 
be  eliminated.  This  consolidation  is  described  in  the  model  by  pre¬ 
scribing  a  consolidation  density  p^  (Point  C  in  Figure  11).  Once  this 


density  is  exceeded,  loading  and  unloading  occur  along  the  solid  P-V 
relation,  curve  ECF  in  Figure  11.  The  solid  curve  in  our  model  is 
described  by  the  relationship 


') 

P  =  Ap  +  Bp- 


(6) 


where 

A  =  bulk  modulus  of  solid  material  at  ambient  pressure 

U  =  p/p  -  1;  p  =  density  of  solid  material  at  ambient  pressure 
so  so 

B  =  coefficient  of  term  describing  monlinear  behavior. 

The  above  pressure-volume  relation  is  consistent  v^;ith  a  linear  relation¬ 
ship  for  bulk  modulus 


if 


K  =  K  +  aP 
o 


a  =  1  + 


2B 

A 


(7) 


The  form  of  Equation  (6)  is  consistent  x>7ith  hydrostatic  pressure-volume 
measurements  made  in  polycrystalline  salt. 

The  material  parameters  required  for  specifying  the  pressure- 
volume  relation  are: 

•  Bulk  modulus  and  density  of  the  porous  material  at 
ambient  pressure  (K^  and  p^) . 

•  Bulk  modulus  and  density  of  the  porous  material  at 

the  consolidation  point  (K  and  p  ). 

c  c 

•  Bulk  modulus  and  density  of  the  solid  material  at 

ambient  pressure  (K  ,  p  )  and  P-V  relation  for  the 
so  so 

solid  material. 

•  ParameLeis  for  the  parabolic  segments  describing  the 
porous  loading  curve  ABC  in  Figure  11. 


The  determination  of  these  parameters  from  the  experimental  data  of 
5  ^ 

Heard  et  al.  is  described  next  and  summarized  in  Table  2. 


Table  2 

PARAMETERS  FOR  POROUS  MATERIAL  MODEL 


Material  parameters 

Bulk  modulus  of  porous  material 

at  ambient  pressure  =  225  kbar  (22.5  GPa) 

Bulk  modulus  at  consolidated 

density  =  274  kbar  (27.4  GPa) 

Coefficient  of  linear  term  in 
the  P-p  relation  for  solid 

material  A  =  245  kbar  (24.5  GPa) 


Parameters  definina  porous  material  mean  stress-volume  path 


Pj^  =  0  or  ambient 
P2  =  0.802  kbar 
Pj  =2.31  kbar 
P^  =5.70  kbar* 
DP^  =  -0.1  kbar 
DP2  =0.12  kbar 
DP^  =  -0.3  kbar 


=  0.4675  cm^/g 
V2  =  0.4652  cm^/g 
=  0.4605  cm^/g 
=  0.4535  cm^/g 


Consolidation  va]ues. 


3 

A  Pj  value  of  2.139  g/cm  was  used  in  our  calculations.  The 
value  was  obtained  from  the  ultrasonic  data  of  Heard  et  al.  (Figure  9 
in  Reference  5).  Beyond  a  confining  pressure  of  25  MPa,  a  straight  line 
can  be  drawn  through  the  K-P  results  obtained  from  their  ultrasonic 


Our  density  and  sound  velocity  measurements  (Section  2)  on  the  salt 

cores  are  similar  to  those  of  Heard  et  al.^  We  will  therefore  assume 
that  all  of  their  measurements  are  applicable  to  our  cores. 
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measurements.  This  straight  line  was  extrapolated  to  P  =  0  and 
provided  a  value  of  =  225  kbar  (22.5  GPa).  At  the  consolidation 
point,  the  bulk  modulus  and  density  of  the  porous  material  are  the 
same  as  for  the  solid  material  and  are  described  below. 


The  initial  bulk  modulus  and  density  of  the  solid  material  are 
the  values  cited  by  Simmons  and  Wang.^^  These  are  =  2A5  kbar  and 


p  =2.16  g/cm"'.  Using  these  value's,  we  fitted  a  quadratic  relation, 
so 

as  given  in  Equation  (6),  to  Decker's  data  on  polycrystalline  salt. 


P  =  245y  +  505. 3y  kbar  (8) 

The  pressure-volume  curve  for  the  porous  solid  is  obtained  by 
fitting  parabolic  segments  to  the  mean  stress-volume  curve  obtained 
under  static  uniaxial  strain  loading.  As  explained  earlier,  this 
procedure  approximates  the  effect  of  shear  deformation  on  compaction. 
Three  parabolic  segments  were  used  to  fit  the  loading  curve;  the  material 
parameters  describing  these  segments  are  presented  in  Table  2.  The 
poro"s  and  solid  loading  curves  used  in  our  work  are  shown  in  Figure  12. 

3 

The  consolidation  point  is  given  by  P^  =  5.7  kbar  and  =  0.4535  cm  /g. 
Beyond  this  point,  the  material  loads  and  unloads  along  the  solid 
curve. 

For  the  porous  material,  the  unloading  and  reloading  path  is 

i* 

specified  by  the  bulk  modulus.  This  requires  values  for  K^.  and 

(bulk  modulus  at  consolidation) .  A  value  of  K  =  274  kbar  is  obtained 

c 

by  using  Equation  (7)  and  the  P^  value  cited  above.  The  bulk  modulus 
at  an  intermediate  point  is  provided  by  linear  interpolation. 


To  get  a  good  match  with  the  porous  data,  we  used  =  2.J.57, 
instead  of  2.16. 

t- 

In  the  porous  region,  the  bulk  modulus  does  not  define  the  loading 
modulus. 
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FIGURE  12  POROUS  AND  SOLID  MATERIAL  MODEL  LOADING  PATH 


DEVIATOR  STRESS  DESCRIPTION 


Heard  et  al.  have  performed  static  uniaxial  stress  and  uniaxial 

5 

strain  experiments  on  salt.  Their  data  provide  the  basis  for  the 
deviator  stress  (or  strength)  model  in  our  work.  We  first  discuss 
their  results  and  then  describe  the  constitutive  model.  A  detailed 
derivation  of  the  equations  and  the  insertion  of  the  deviator  stress 
subroutine  into  the  COPS  code  is  presented  in  Appendix  C. 

Data  from  uniaxial  stress  experiments  show  a  brittle-ductile 
transition  at  a  confining  pressure  of  10  bars.^  Beyond  this  value, 
the  increase  in  strength  with  pressure  is  negligible.  During  loading, 
the  material  displays  yield  followed  by  strain  hardening.  Unloading 
in  these  experiments  is  elastic  and  results  in  a  small  amount  of  net 
compaction  caused  by  the  small  porosity  present  in  the  salt.  Unlike 
other  geologic  solids  of  interest,  salt  deforms  like  a  metal  with  the 
following  characteristic  features:  distinct  yield,  nonlinear  strain 
hardening,  negligible  dependence  of  yield  strength  of  pressure,  and 
no  dilatancy. 

The  uniaxial  strain  data  plotted  by  Heard  et  al.^  show  increasing 
shear  strength  with  confining  pressure  when  plotted  in  stress  space. 
This  strength  increase  is  deceptive  because  the  deviatoric  strains 
are  also  increasing  with  compression.  When  transformed  to  a  deviator 
stress-strain  plot,  these  data  show  strain  hardening  somewhat  similar 
to  that  observed  in  the  uniaxial  stress  experiments. 

Based  on  these  result.,,  the  following  elastic-plastic  model  with 
isotropic  strain  hardening  is  used  to  represent  the  deviator  stress- 
strain  relations. 

(1)  Additivity  of  incremental  elastic  and  plastic  strains 


de. .  =  de.®  +  de.P 

ij 


(9) 


This  classification  also  includes  tests  in  which  the  axial  stress  is 
superposed  on  a  constant  confining  hydrostatic  pressure. 


(2)  Incremental  Hooke's  law 


do. .  “  2Gdc , . 
ij  ij 


where  de . T 
ij 


de. .  - 


(3)  Yield  surface  of  the  von-Mises  type  with  strain  hardening 


f  =  v/j'  -  Y  =  0 


where  =\  ^  O.T  •  a.' 

V  2  \/  2  ij  ij 


and  Y  =  Y^  +  g(e?j) 


(4)  Plastic  incompressibility 


de'^  =  0 


G  is  the  shear  modulus  and  is  a  function  of  pressure;  g  is  a  strain¬ 
hardening  function  to  be  determined  from  the  experimental  data. 

The  above  equations  are  the  usual  relations  presented  in  a 
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continuum  plasticity  textbook.  Given  the  available  experimental 
data  on  salt,  these  equations  are  quite  adequate  as  will  be  seen  in  the 
next  section. 

The  above  equations  represent  a  model  that  is  suitable  for 

describing  static  tests.  For  simulating  impact  experiments,  we  need 

to  incorporate  a  rate  dependent  or  viscoplastic  model.  The  model  used 

13 

here  is  similar  to  the  overstress  model  discussed  by  Perzyna.  The 
stresses  beyond  yield  are  given  by 


o'..  =  (o:.)^  1  - 

ij  13  I 


The  subscripts  I  and  E  denote  the  instantaneous  and  equilibrium  values 
as  discussed  below;  At  and  are  the  time  step  and  relaxation  time, 

also  explained  in  the  following.  In  a  rate- independent  problem,  stress 
states  outside  the  yield  surface,  that  is,  f  >  0,  are  not  permissible. 
The  magnitude  of  the  plastic  strain  is  determined  by  the  condition 
that  df  =  0.  For  a  viscoplastic  model,  stress  states  f  >  0  are 
permissible  and  the  stresses  relax  torward  the  current  value  of  the 
static  yield  surface  Y  (equilibrium  value)  with  a  relaxation  function 


proportional  to  the  stress  overshoot  and  a  relaxation  time  T^-  The 
instantaneous  stress  is  defined  as 


(a..)'  =  (a. .)'  +  2G  Ae.: 
'  xj  I  ij'o  IJ 


where  (a  '.)  is  the  deviator  stress  at  the  start  of  the  increment  and 
°  * 

Ae.'  is  the  deviatoric  strain  increment  in  a  time  step  At.  The 
corrected  stresses  at  the  end  of  the  increment  are  given  by  0^7,  and 
the  plastic  strain  increments  are  given  by 

(0.:),  -  o.: 

For  At/T^  «  1,  Equation  (13)  gives  0^'  s  ^t/T^  =  1, 

0.'  =  (0.r)_.  For  At  >  T„,  the  stresses  are  set  equal  to  (0.O„. 
ij  ij  E  r’  ^  xj  E 

Related  forms  for  Equation  (13)  may  be  seen  in  studies  by  Seaman^^  and 

Johnson. 

For  a  perfectly  plastic  material  (including  pressure-dependent 
strength),  the  use  of  Equation  (13)  is  quite  straightfon^ard  because 
the  yield  surface  is  fixed  in  stress  space.  In  the  presence  of  strain 
hardening,  the  determination  of  0^^  and  is  complicated  because  of 

changing  Y.  However,  a  closed  form  solution  within  an  increment  has 
been  obtained.  Details  of  this  derivation  are  given  in  Appendix  C; 
here  v^e  present  the  results  of  this  derivation. 


Assumed  elastic  for  computing 
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a> 

I  I 


A 


i 

fm 


We  have  chosen  a  strain-hardening  model  given  by 


Y  =  Y  +  M  VAi; 


where  >/  AI^  ^  =  ^/l/2  cT?  •  AeTP  and  Y  is  the  yeild  stress  at  the 

2  ij  ij  o 

start  of  the  increment-  M  is  a  nonlinear  function  and  depends  on  the 

current  value  of  the  accumulated  plastic  strain  e.P.  in 

ij  2E 

Equation  (13)  is  the  same  as  Y. 

The  first  terra  on  the  right  side  is  the  equilibrium  stress  at  the 
start  of  the  increment.  The  above  equations  can  be  combined  to  give 


2  2GT 

M  +  ■  ^ 

At 


Because  all  the  terras  on  the  right  side  of  Equation  (18)  are  known, 
we  can  obtain  Then  Equations  (17),  (13),  and  (15)  can  be 

evaluated  to  give  all  the  stresses  and  strains. 

The  material  parameters  for  the  deviator  stress  model  are: 

•  Shear  modulus  G  and  its  dependence  on  pressure  P. 

•  Initial  value  of  yield  stress  and  the  strain  hardening 
relation. 

•  Relaxation  time,  T  . 

r 

We  consider  the  first  two  items  here  and  postpone  determination  of 

T  to  Section  *. 
r 

The  shear  modulus  was  assumed  to  vary  linearly  with  pressure. 


G  =  G  +  G-  P 
o  1 


For  G  ,  we  used  a  value  of  13.7  GPa  (137  kbar)  given  by  Heard  et  al. 
(Table  3  of  Reference  5).  A  Gj^  value  of  0.22  GPa  (2.2  kbar)  was 


. . . ililii . MW . I! . IB . . . . . . . . . . . . . . . . . . . 


determined  by  taking  the  slope  of  the  shear  modulus-pressure  curve 
(Figure  9  in  Reference  5) . 


An  initial  yield  stress  value,  =  25.4  MPa,  was  obtained  from 
the  quasi-static  uniaxial  stress  data  (Figure  4  in  Reference  5). 

Stress  difference  (a^  -  o^)  versus  axial  strain  (e^^)  and  radial  and 
hoop  strains  (e^  -  e^)  curves  obtained  from  these  data  were  used  to 
construct  the  strain-hardening  curve  (/J^  -  Y^  versus  shown  in 

Figure  13.  This  curve  was  fitted  by  a  series  of  parabolic  segments 
in  a  manner  similar  to  that  used  for  the  porous  model.  The  parameters 
corresponding  to  the  six  segments  are  given  in  Table  3.  Segment  6  and 
part  of  segment  5  are  not  shown  in  Figure  13  because  the  strain¬ 
hardening  curve  in  this  region  is  nearly  a  straight  line. 

Table  3 

PARAMETERS  DEFINING  THE  STRAIN-HARDENING  CURVE 


/j^  -  Y^  (kbar) 

A[/J2e  ■  (kbar) 

0.0 

0.0 

0.0 

0.05197 

0.007 

0.000965 

0.10393 

0.01366 

0.0042395 

0.195 

0.004 

0.021223 

0.218 

0.0 

0.037882 

0.22 

0.0 

0.068 

0.22 

0.00 

1.0 

Al.^  -  Y^]  is  the  maximum  deviation  from  a 
zt.  o 

linear  /Jo_  -  Y  versus  /I^^  variation. 

2E  o  2 


SIMULATION  OF  STATIC  DATA 

The  strength  model  parameters  presented  in  the  previous  subsection 
were  determined  from  the  uniaxial  stress  data.  Hence  no  attempt  was 
made  to  simulate  the  uniaxial  stress  experiments.  Instead,  we  simulated 
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FIGURE  13  STRAIN  HARDENING  CURVE  (SOLID  LINE)  OBTAINED  FROM  THE  DATA 
OF  HEARD  ET  AL. 

The  dotted  line  shown  in  segments  3  and  5  is  the  curve  actually  used. 


the  uniaxial  strain  data.  This  simulation  provides  a  good  check  of  the 
strength  model  because  it  exercises  the  model  for  a  different  test 
condition. 

The  uniaxial  strain  experiment  is  simulated  using  a  test  program 
that  imposes  one-dimensional  strain  increments  to  the  desired  peak 
strain  value.  From  a  knowledge  of  the  strains,  the  stresses  are 
calculated  using  the  material  model  routines  POROUSl  (P-V  relation) 
and  VPSH  (deviator  stress-strain  relation) .  The  listings  for  these 
routines  are  given  in  Appendices  B  and  C.  This  simulation  is  similar 
to  the  use  of  material  property  subroutines  in  the  wave  propagation 
code. 

Figure  lA  compares  the  calculations  with  the  experimental  measure¬ 
ments  of  Heard  et  al.^  Figure  lA(a)  shows  the  plots  for  versus 
V/V^,  and  Figure  lA(b)  shows  the  plots  for  versus  o^.  The 

agreement  is  generally  quite  good  except  over  the  stress  range  in  which 
pore-collapse  predominates.  The  experimental  data  show  that  a  strength 
decrease  accompanies  pore-compaction.  However,  in  our  model  no  such 
interaction  exists  between  pore-compaction  and  material  strength. 
Because  of  the  complexity  in  describing  this  interaction,  no  attempts 
were  made  to  include  it  in  the  material  model.  As  remarked  earlier, 
the  increase  in  with  increasing  is  misleading  because  the 

real  reason  for  this  increase  is  strain  hardening  and  not  pressure 
dependence.  From  this  simulation,  we  conclude  that  the  strength  model 
is  satisfactory  for  calculations  along  other  load  paths. 


MPa 


0.960  0.965  0.970  0.975  0.980  0.985  0.990  0.995  1 

V/Vq 

(a)  o,  VERSUS  V/Vq 


FIGURE  14  COMPARISON  OF  THE  CALCULATED  AND  EXPERIMENTALLY 
MEASURED  CURVES  FOR  STATIC  UNIAXIAL  STRAIN  LOADING 
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ANALYSES  AND  CONCLUSIONS 


4 

-  -9 


In  this  section  ».'e  analyze  the  data  presented  in  Section  2.  First 
the  Jongitudinal  profiles  are  analyzed  to  give  the  longitudinal  stress- 
volume  results.  This  is  followed  by  a  numerical  simulation  of  the  wave 
profiles  and  comparison  with  experimental  data.  Next,  the  rapid  decay 
of  shear  wave  amplitude  is  discussed  by  analyzing  the  equations  governing 
wave  propagation.  Some  concluding  remarks,  including  directions  for 
future  work,  are  presented  at  the  end  of  the  section. 

LONGITUDINAL  STRESS-VOLUME  RELATION 

In  a  previous  report  we  described  the  devalopuseut  of  a  Lagrangian 
analysis  for  combined  compression  and  shear  waves.  The  analysis  consists 
of  solving  the  governing  equations  (Equations  A. 5  through  .A.S  in  Appendix 
A)  from  a  knowledge  of  U2^(t)  and  u,(c)  profiles  at  several  locations. 

It  has  been  shown  previously  that  the  longitudinal  stress-volume  (a^,-P) 
relation  is  completely  determined  by  the  Uj(t)  profiles.  This  is  true 
whether  we  have  biwaves  or  uniaxial  strain  waves.  The  main  advantage 
of  the  Lagrangian  analysis  is  that  the  data  are  cast  in  a  form  th.at  is 
easier  to  use  for  developing  constitutive  nsodels.  The  results  of  an 
analysis  of  the  u^f.t)  profiles  are  given  below. 

The  uniaxial  strain  profiles  from  experiment  4  (78-2-25)  vrere 
analyzed  to  give  the  stress-volume  loading  and  unloading  path  shown 
in  Figure  15.  .An  elastic  loading  to  1  kbar  can  be  observed.  .A  peak 
stress  of  10.6  kbar  indicates  the  Hugoniot  state.  The  loading  and 
unloading  states,  except  for  late  time  unloading,  are  nearly  similar 
at  all  gage  locations. 

A  similar  analysis  was  performed  for  experiment  3  (78-2-27)  to 
provide  the  longitudinal  stress-volume  data  under  compress* on  and  shear 
loading.  The  loading  path  and  peak  states  from  experiments  3  and  4 
are  plotted  with  the  LLL  ds ta  in  Figure  16.  The  broken  line  is  the 
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LLL  static  uniaxial  strain  data,’  and  the  circle  points  are  the  LLL 
Hugoniot  data.  The  LLL  Hugoniot  data  indicate  a  more  compliant  material 
response  chan  our  data.  We  believe  that  this  discrepancy  is  caused  by 
differences  in  the  material  in  our  work  and  the  LLL  dynamic  work.  The 
density  and  longitudinal  sound  speed,  reported  by  Anderson  and  Larson,^ 
are  smaller  by  0.5%  and  5%,  respectively.  Although  all  the  differences 
are  close  to  experimental  accuracy  limits,  the  consistently  lower 
values  for  each  case  indicate  that  these  differences  are  real.  Our 
density  and  sound  speed  measurements  are  close  to  the  LLL  static  work. 

The  main  result  inferred  from  Figure  16  is  that  the  initial  loading 
modulus  is  higher  for  the  dynamic  data  than  for  the  static  data.  At 
higher  compressions,  however,  an  extrapolation  of  the  static  data 
agrees  with  the  final  states  observed  in  our  work.  The  higher  loading 
modulus  at  lower  stresses  is  a  consequence  of  rate-dependent  pore- 
collapse;  that  is,  for  a  given  stress,  less  pore-compaction  is  observed 
at  high  rates  than  at  low  rates. 


H  * 


NUMERICAL  CALCULATIONS  OF  WAVE  PROFILES 

To  assess  the  constitutive  model  for  salt  presented  in  Section  3 

and  to  extend  it  for  handling  dynamic  material  response,  we  made  numerical 

simulations  of  the  observed  wave  profiles,  using  the  wave  propagation 
'  16 

program,  COPS.  ’  The  COPS  code  was  developed  to  study  one-dimensional 
compression  and  shear  wave  propagation  in  solids.  It  is  an  extension 

14 

of  the  one-dimensional  codes  used  for  studying  uniaxial  strain  waves. 

The  impacting  solids  are  modeled  as  planar  slabs,  and  the  governing 
equations  are  solved  using  the  artificial  viscosity  method  pioneered 
by  von  Numann  and  Richtmeyer.^^  To  handle  shear  wave  propagation, 
the  code  must  also:  (1)  incorporate  material  rotation  and  consistent 
use  of  stress  and  strain  measures,  (2)  use  artificial  viscosity  for 
shear  waves,  (3)  minimize  the  dispersion  in  the  shear  wave  caused  by 
the  artificial  viscosity  and  (4)  properly  specify  the  shear  particle 
velocity  at  the  impact  surface. 

The  material  model  parameters,  derived  from  the  static  parameters, 
were  presented  in  Section  3.  Adaptation  of  the  model  to  the  dynamic 
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calculations  requires  the  specification  of  a  relaxation  time  in  the 
deviator  stress  model  as  indicated  in  Section  3  and  a  slightly  stiffer 
P-V  relation.  Both  of  these  changes  were  incorporated  by  simulating 
the  LLL  uniaxial  strain  experiment  having  a  peak  particle  velocity 
of  0.138  mm/ps  (experiment  NN-6  in  reference  6).  The  values  for  the 
relaxation  time  and  the  dynamic  parameters  for  the  P-V  relation 
are  listed  in  Table  4.  All  the  simulations  presented  below  were  made 
with  the  same  parameters. 


Table  4 

DYNAMIC  PARAMETERS 


Pj^  =  0  kbar 

A 

P2  =  0.5  kbar 

A 

P^  =  2.8  kbar 

P,  =5.7  kbar 
4 

DPj^  =  0  kbar 
DP2  =  0.12  kbar 
DP^  =  -0.3  kbar 


0.05  ysec 

=  0.4675  cm^/g 
V2  =  0.4663  cm^/g 
Vj  =  0.4605  cm^/g 
=  0.4535  cm^/g 


These  are  the  values  that  are  different 
from  the  corresponding  values  shown  in 
Table  3. 


Figure  17(a)  compares  the  numerically  calculated  wave  profiles 
with  the  experimental  data  from  LLL  experiment  NN-6.  The  agreement  is 
quite  good  except  that  the  calculated  plastic  wave  is  slightly  stiffer. 
The  irregular  profile  in  gage  2  is  an  experimental  artifact  (either  in 


A 

In  our  work  the  dynamic  relation  was  chosen  by  matching  the  wave 
arrival  in  experiment  NN-6  of  reference  6. 
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measurement  or  in  data  reduction)  and  was  not  observed  in  any  other 
record. 

Figure  17(b)  shows  the  comparison  for  a  higher  stress  experiment 
(NN-7  in  reference  6).  Good  agreement  is  again  observed  between  the 
experimental  data  and  calculations.  The  slight  increase  in  the  particle 
velocity  data  at  late  time  is  an  experimental  artifact  caused  by  lead 
stretching.  In  this  comparison,  the  experimental  profiles  are  steeper. 

These  two  simulations  show  that  the  model  presented  here  adequately 
fits  the  shock  data  obtained  on  salt.  Using  the  same  model,  we  simulated 
the  combined  compression  and  shear  wave  experiments  in  a  similar  manner. 

,  Particle  velocities  measured  at  the  impact  surface  during  experiments 

A 

3  and  4  (Section  2)  were  prescribed  and  wave  profiles  were  obtained  at 
!  distances  comparable  to  the  gage  locations.  To  minimize  shear  wave 

dispersion  caused  by  artificial  viscosity,  we  selected  the  lowest- 
allowable  value  of  artificial  viscosity. 


Figure  18(a)  compares  the  experimentally  measured  longitudinal 
profiles  with  the  calculated  wave  profiles.  Except  for  the  steepness 
near  the  top  of  the  plastic  wave,  the  agreement  is  quite  good.  The 


numerically  calculated  and  experimentally  measured  shear  wave  profiles 
are  shown  in  Figure  18(b).  The  calculated  wave  profiles  show  wave 
arrival  with  an  elastic  shear  wave  speed  followed  by  a  very  gradual 
increase  in  particle  velocity.  The  peak  amplitudes  are  much  smaller 
than  the  prescribed  particle  velocity  amplitude;  that  is,  the  shear 
wave  amplitude  is  rapidly  attenuated  near  the  impact  surface.  This 
result  is  qualitatively  similar  to  experimental  observations.  However, 
the  arrival  time  and  the  structure  in  the  shear  wave  data  are  very 
different  from  the  calculated  profiles.  It  is  quite  certain  that  the 


faster  traveling  longitudinal  wave  produces  a  voltage  signal  because 
of  experimental  artifacts  (Section  2),  but  the  exact  magnitude  of  this 
signal  is  not  easily  known  and  cannot  be  easily  subtracted  from  the 


In  these  simulations  we  need  to  simultaneously  prescribe  Uj^(t)  and 
u^(t)  at  the  impact  surface. 
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observed  data.  If,  in  Figure  18(b),  we  subtract  the  particle  velocity 

* 

amplitude  in  the  experimental  data  at  times  corresponding  to  U.8  ys 
and  1.8  ys  from  the  late  time  amplitude,  the  remaining  amplitude  values 
are  comparable  to  the  calculated  amplitudes.  Unfortunately,  there  is 
little  justification  for  this  procedure. 

The  compression  and  shear  wave  simulations  show  that  the  consti¬ 
tutive  model  predicts  the  longitudinal  wave  profiles  reasonably  well. 

The  calculated  shear  wave  profiles  show  rapid  attenuation  with  propa¬ 
gation  distance,  and  the  calculated  amplitudes  are  small  in  accordance 
with  experimental  data.  However,  the  arrival  time  and  structure  of  the 
observed  profiles  are  quite  different  from  the  calculated  profiles. 

The  probable  causes  of  this  discrepancy  are  experimental  artifacts;  it 
is  not  possible  to  use  the  shear  nave  data  for  an  improvement  of  the 
material  model  outlined  in  Section  3. 


LOCALIZATION  OF  THE  SHEAR  WAVE 

The  shear  wave  data  on  salt  showed  extremely  rapid  attenuation 
near  the  impact  surface  and  very  little  shear  motion  coupled  to  the 
specimen  interior.  A  similar  result  was  also  obtained  from  the  numerical 
calculations.  This  rapid  shear  wave  attenuation  with  propagation 
distance  can  be  explained  by  a  simple  wave  analysis  as  shown  below. 

The  equations  describing  shear  x^ave  propagation  can  be  written 
as  (see  Appendix  A  and  reference  4) 


ay  = 

9t  2  9  h 


(21) 


These  are  times  corresponding  to  the  arrival  of  the  shear  waves  in 
the  calculations. 


where  we  have  written  h  =  Y  =  £^^2  ^21^’  ^  ~  ^12  convenience 

of  typing.  We  now  define  a  velocity 


:  -  (— '1 


Inhere  the  symbol  a  denotes  either  T,  Y»  or  U2.  Hence,  represents  a 
velocity  associated  with  a  particular  a  value  of  the  wave.  Definition 
(22)  is  identical  to  the  contour  velocities  introduced  by  Fowles  in  the 
analysis  of  uniaxial  strain  leaves. 

Recognizing  that  U2  is  a  function  of  the  independent  variables  h 
and  t  we  can  derive 


C  = 

“2 


-  Ou2/3t)^ 
('dU2/3h)^ 


CcTbining  Equations  (21)  and  (23),  we  can  express  the  change  in  y  at 
a  gage  location  as 


dY  =  — along  h  =  constant 


Knowledge  of  C  as  a  function  of  u.^  enables  us  to  integrate  Equation 

U2 

(24)  and  determine  Y(t)*  An  interesting  aspect  of  Equation  (24) 
pertinent  to  the  present  work  concerns  the  relative  magnitudes  of  du, 
and  C  .  For  du~  »  2C  ,  a  very  large  increment  in  strain,  dY,  is 

”2  ^  ”2 

obtained.  Unlike  the  compiession  wave,  this  inequality  is  very  plausible 
for  the  shear  wave.  The  consequences  of  this  result  are  explored  next. 


u.,  =  (h,c) 


<'“2  -  {-^)  ^  (•^) 

For  constant  U2  we  get  equation  (23), 
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Defining  similar  to  and  combining  with  Equations  (20), 

(21) ,  and  (23) ,  can  write 

dT 

— —  =  2p  C  C  along  h  =  constant  (25) 

dy  o  T  ^ 

Equation  (25)  expresses  the  modulus  or  the  slope  of  the  shear  stress- 
strain  relation  at  a  gage  location.  Knowledge  of  and  over  the 

h-t  domain  of  an  experiment  can  be  used  to  obtain  the  dynamic  T-y 
relation.  A  different  but  related  procedure  is  used  in  the  Lagrangian 
analysis  presented  in  reference  A. 

Equations  (2A)  and  (25)  are  applicable  in  general.  To  derive 

the  relationship  we  seek,  we  will  make  some  simplifying  assumptions 

about  the  material.  If  we  assume  chat  T=  T(y)  and  that  wave  propagation 

18 

is  into  a  constant  region,  as  is  generally  the  case,  we  can  write 


(26) 


Equation  (26)  can  be  combined  with  Equations  (25)  and  (2A)  to  give 


(l7)h 


and 


du. 


dy  =  - 


(27) 


(28) 


Equation  (28)  shows  that  dy  is  determined  by  the  relative  magnitudes 

of  du2  and  the  shear  stress-strain  modulus.  That  is,  large  strains  can 

build  up  (also  termed  ’localization’)  on  the  microsecond  time  scale  even 

if  C  is  significantly  positive.  Because  this  localization  is  controlled 
“2 

in  part  by  the  imposed  motion,  we  refer  to  it  as  ’kinematic  Ir oaliza'^ion’ . 


Unlike  compressive  deformation,  large  shear  strains  inhibit 
propagation  of  large  amplitude  shear  waves.  This  happens  because  large 


shear  strain  results  in  rupture  (resulting  in  stress  drop)  or  a 
flattening  of  the  shear  stress-strain  curve.  The  second  effect  occurs 
in  our  calculations,  with  the  shear  wave  being  localized  near  the 
impact  surface.  A  similar  localization  is  generally  not  observed  xn 
compressive  wave  propagation  because  the  loading  curve  is  concave 

upwards  due  to  the  curvature  in  the  pressure-volume  relation. 

CONCLUSIONS 

The  objectives  of  the  present  work  were  to  obtain  the  compression 
and  shear  wave  profiles  in  salt  and  to  interpret  these  data  in  terras  of 
a  constitutive  moael.  These  objectives  have  largely  been  met  except  for 
determining  th2  precise  cause  of  the  observed  structuxe  in  the  shea* 
wave  profiles.  With  further  effort,  these  causes  could  probably  be 
determined  and  mitigated.  The  benefits  of  such  an  effort,  from  a 
practical  viewpoint,  are  not  obvious.  The  reasonable  agreement  between 
the  simple  constitutive  model  developed  here  and  the  experimental 
profiles  suggests  that  similar  moat.'’ing  on  natural  salt  is  warranted. 
Furthermore,  emphasis  should  be  placed  on  accurate  determination  and 
modeling  of  the  unloading  portion  of  the  longtitudinal  profiles. 

The  granite  results,  though  from  a  single  experiment,  suggest 
several  areas  of  investigation  in  hard  rocks.  It  would  be  interesting 
to  determine  if  the  amplitude  of  the  shear  wave  that  can  be  propagated 
into  the  material  is  related  to  the  amplitude  of  the  compression  wave. 
Tne  answer  to  this  question  depends  on  two  competing  effects  and  is 
not  obvious.  Increasing  increases  both  the  mean  stress  and  stress 
deviators.  The  increase  in  mean  stress  favors  a  larger  amplitude  shear 
wave,  v’hereas  an  increase  in  stress  deviator  favors  a  reduction  in  the 
shear  amplitude. 

Because  the  shear  wave  propagates  into  a  nearly  constant  region 

19 

in  the  granite,  the  method  developed  by  Gupta  et  al.  can  be  used  to 

determine  the  dynamic  mean  stress-volume  (a  / 3  versus  V/V  )  relation. 

mm  o 

In  this  method  a  bulk  modulus-density  (K-p)  relation  is  determined  from 
experimental  measurements  of  longitudinal  and  shear  wave  speeds  in  the 


shocked  state.  The  K-p  curve  is  integrated  to  obtain  the  a  -p/p 

mm  o 

relation.  The  knowledge  of  mean  stress-volume  relation  permits  a 

determination  of  the  dynamic  strength  from  existing  uniaxial  strain 

data.  It  would  be  desirable  to  determine  if  the  material  loses  strength 

(under  uniaxial  strain)  beyond  a  threshold  compressive  stress.  In 

PMMA  the  dynamic  material  strength  decreased  beyond  a  threshold  stress 

giving  opposite  results  under  dynamic  and  static  uniaxial  strain 
19 

ixperimenls. 

Finally,  wave  propagation  across  faults  and  joints  is  of  interest 
to  underground  testing.  An  important  question  in  these  studies  concerns 
the  5>agnitude  of  dynamic  friction  and  post-sliding  behavior,  labile 
phene menology  experiments  at  small  scale  are  useful,  the  inclined  impact 
experiments  can  provide  th*'  quantitative  data  needed  for  material 
modeling.  It  would  be  usefi’l  to  quantify  the  frictional  behavior  across 
interfaces  in  the  presence  of  inelastic  deformation  because  present 
knowledge  in  this  area  is  quite  limited.  A  direct  measure  of  the 
frictional  behavior  at  an  interface  is  determined  from  the  shear  particle 
velocities  across  an  interface. 
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Appendix  A 

COMPRESSION*  AND  SHEAR  WAVES  IN  IMPACTED  SOLIDS 

dfc 

This  appendix  sucsaarizes  the  IMPS  sethod  for  coapression  and  shear 
vave  rseasurenents  and  presents  the  appropriate  governing  equations  for 
describing  wave  propagation.  An  error  analysis  for  shear  particle 
velocity  seasureisents  is  also  presented. 

EXPERIME-NTAL  ME.ASUREMENTS 

The  experinentai  aethod  is  based  on  the  neasuresent  of  particle 
velocities,  using  Faraday’s  lav  for  coving  circuits,  in  the  inclined 
plate  icpact  experisent  shown  in  Figure  .A.I.  .A  flyer  plate  with  velocity 
along  the  X'-axis  icpacts  a  parallel  target  plate.  The  norcal  to  the 
icpacting  surfaces  is  along  the  Xj^-axis.  Upon  icpact,  the  target  plate 
is  iaparted  a  cot ion  having  coaponents  along  the  X^  and  X,  axes.  By 
inserting  chin  csetallic  gages  in  the  X^-X^  plane  and  censuring  the  EMF, 

£,  generated  by  their  notion  in  a  constant  cagnetic  field,  we  can  deteraine 


c’ne  particle  velocities  free  the  relation 


<u  X  B) 


(A.l) 


where  t  is  the  length  vector  for  the  gage,  u  is  particle  velo«  ey,  and 
B  is  the  cagnetic  field.  By  ceasuring  E  for  two  different  B  field 
orientations,"  we  can  obtain  the  particle  velocity  cise  profiles  UjCt) 
and  u^(c).  Icproved  .accuracy  is  obtained  fay  ceasuring  either  u, (t)  or 

i  * 

u,(t)  profiles  in  a  single  experiaent. * 


J_nCemal  Measarecents  of  and  S  Waves. 

‘E:-;cepc  for  the  B  fields,  the  e.eperisents  should  be  id^-ntical. 

' 1 t  has  been  shown  elsewhere  that  expericents  can  be  designed  so  that 
Uj(c)  and  u,(t)  can  be  ceasured  sicultaneously  in  a  single  experisent. 

However,  this  has  not  been  done  to  date. 


a  II)  i*<  . . .  . . .  . . .  i  ^ 


An  experimental  facility  for  impacting  parallel  inclined  plates 

•k 

was  designed  and  constructed  under  a  previous  DNA  contract.  A  magnet 
system  was  developed  that  provides  the  desired  orientation  of  the  B 
field  in  the  X2^-X2  plane.  Target  construction  for  these  experiments 
is  similar  to  the  usual  uniaxial  strain  experiment.  In  addition  to  the 
particle  velocity  profiles,  each  experiment  measures  impact  velocity, 
tilt,  and  the  magnetic  field  strengths.  The  errors  in  shear  particle 
velocity  measurements  are  discussed  below. 

Because  the  amplitude  of  the  shear  signal  in  the  salt  experiments 
was  very  small,  it  was  important  to  establish  the  accuracy  of  this 
signal.  To  achieve  this  objective,  we  examined  the  magnitude  of  the 
signal  that  may  be  caused  by  the  preceding  longitudinal  wave.  Note 
that  we  are  not  considering  any  material  effects  such  as  coupling  of 
longitudinal  and  shear  motions  due  to  nonlinear  material  response; 
instead,  we  are  concerned  only  with  experimental  artifacts. 

The  particle  velocity  for  the  longitudinal  wave  is  given  by  u  = 

(u^,  0,  0).  For  B  =  (Bj^,  0,  0)  and  £  =  ^3)*^  the  longitudinal 

wave  does  not  lead  to  a  signal  under  ideal  conditions  (as  determined 
from  Equation  A.l).  However,  all  the  above  quantities  can  have  small 
errors  because  of  impact  tilt,  field  misalignment,  and  gage  placement. 

We  now  estimate  the  error  in  the  signal,  Ae,  caused  by  experimental 

errors  Au,  AB,  and  A£  associated  with  the  longitudinal  wave.  We  will 

7 

ignore  terms  of  the  order  of  A"  because  A  itself  is  a  very  small  quantity. 

-h  - 

We  denote  u  and  B  as  follows:  u  =  (u^,  U2,  and  B  =  (B^,  B2>  ^3/* 

Here  Au.  (or  AB.)  denotes  the  ith  component  of  the  small  error  in  u  or 
B.  We  can  then  write 

Ae  =  ~  ®1  ^^2^  (A. 2) 


Y.  M.  Gupta,  "Development  of  a  Method  for  Determining  Dynamic  Shear 
Properties,"  SRI  International,  Menlo  Park,  CA,  Draft  Final  Report 
Submitted  to  the  Defense  Nuclear  Agency  (1978). 

4.  >> 

‘The  use  of  £  =  (£  ,  £2*  considers  the  entire  gage  and  not  just 
the  active  element. 
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Figure  A. 2  shows  the  gage  layout.  To  obtain  we  must  evaluate 

Equation  (A. 2)  for  each  of  the  three  segments  AB,  BC,  and  CD.  Each 

of  these  segments  is  supposed  to  be  parallel  to  either  the  X2  or 

axis,  but  it  may  have  small  deviations  because  of  incorrect  placement. 

However,  the  terms  arising  because  of  the  incorrect  placement  are  again 
2 

of  the  order  of  A  and  can  be  neglected.  We  can  therefore  write 


In  Equation  (A. 3),  and  have  opposite  signs.  Thus  the  first 
term  on  the  right  side  of  Equation  (A. 3)  is  again  of  the  order  of 
and  is  neglected.  Hence, 


(A. 4) 


Because  the  errors  in  u  and  B  are  independent,  we  need  to  add  the  two 
terms  in  the  parenthesis.  The  error  AB2  is  measured  in  each  experiment. 
The  error  Au2  arises  because  of  tilt  and  is  estimated  by  determining 
the  inclination  of  the  wave  front  from  tilt  records. 

Some  typical  values  pertinent  to  Experiment  2  (78-2-33)  described 
in  Section  2  are 


Jlgl,  =  7.5  mm 
Uj^  =0.2  mm/ys 
Bj^  =  3600  gauss 

AB2  =  50  gauss  (for  this  particular  shot) 

_3 

Au.,  =  1.4  X  10  mm/ys  (taking  the  worst  case) 


Using  the  above  values  we  obtain  Ae,_  =  11  mV.  This  value  of  Ae.„  is 

AD  AD 


equivalent  to  a  Au2  value  given  by 


Au2  =  2  X  10  Uj^ 


(A.5) 
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With  more  care,  these  errors  can  be  reduced  by  a  factor  of  2  or  3. 
However,  for  the  present  experiments  the  value  of  U2  in  Equation  (A. 5) 
represents  the  accuracy  limit  of  the  shear  wave  amplitude  in  a  combined 
compression  and  shear  wave  experiment. 

GOVERNING  EQUATIONS 

In  describing  the  wave  propagation  for  the  impact  situation  shown 
in  Figure  A.l,  it  is  convenient  to  use  the  X^-system.  The  governing 
equations  are  then  one-dimensional;  that  is,  variations  with  respect 
to  only  need  to  be  considered.  The  current  configuration,  particle 
velocity,  and  deformation  gradient  are  defined,  respectively,  as 


X.  =  x.(X  ,t) 

1  1  m’ 

(A.6) 

“1  '(  \ 

m 

(A.7) 

'» -tt) 

(A.  8) 

With  these  definitions, 
propagation  as 


we  can  write 

the  equations 

fr) 

••ra  ■ 

(5) 

(-":■)  ■ 

governing  wave 

(A.9) 

(A.IO) 

(A. 11) 


(A.12) 


where 


T, .  =  First  Piola-Kirchoff  stress  tensor 
li 


=  initial  density 
V  =  specific  volume 


E.  .  = 
ij 


,  /  9x  3x 
1  / _ 

2  \  9x.  ax. 

\  1  J 


-  6..  1  ,  a  finite  strain  measure 


specific  internal  energy. 


For  the  one-dimensional  problem  under  consideration  here,  the  stresses 
and  are  equal  to  the  Cauchy  stresses  and  While  0^2 

is  equal  to  O21’  ^Q^al.  The  use  of  the  symbol 

T  in  Equation  (A. 9)  emphasizes  this  difference. 


Appendix  B 


SUBROUTiME  POROUS! —POROUS  P-V  MODEL 


POROUSl  acts  as  a  porous  material  P-V  relationship-  It  consists 
of  two  sections:  one  section  reads  in  the  required  data  and  initializes 
necessar>'  parameters,  and  a  second  section  computes  the  pressure 
corresponding  to  density  increments  determined  by  the  wave  propagation 
code.  POROUSl  can  be  inserted  into  one-,  two-,  'r  three-dimensional 
codes. 

Table  B.l  is  a  glossary  of  the  main  variables  used  in  POROUSl. 

A  flow  chart  in  Table  B.2  outlines  how  the  equations  presented  in 
Section  3  are  used  in  the  numerical  calculations. 


The  argument  list  of  POROUSl  is 

SUBROUTINE  POROUSl  (NCALL,  IN,  JCELL,  M  M4TL,  DI,  EMGl, 

EMG2,  EJ!G3,  EMG4,  D,  DOLD  EINT,  EOLD, 

P,  NCYC). 

Most  of  these  parameters  are  defined  in  Table  B.l. 

NCALL  directs  control  to  the  data  input  and  parameter  initialization 
section  or  to  the  pressure  computation  section  of  POROUSl.  POROUSl 
data  input  and  parameter  initialization  are  performed  as  part  of  the 
problem-generating  portion  of  the  code.  For  example,  SUBROUTINE  SETUP 
is  the  primary  problem-generating  portion  of  the  COPS  code.  Before 
the  CALL  to  POROUSl,  NCALL  and  IN  are  set  e.u  al  to  0  and  5,  r''.spectively, 
indicating  that  control  is  to  be  transferred  to  the  data  input  and 
parameter  initialization  section  of  POROUSl  and  that  the  data  are  to 
be  read  in  from  input  file  5.  The  CALL  to  POROUSl  in  SETUP  is 

CALL  POROUSl  (NC.4LL,  /.N,  .IE,  M,  MATL,  RHO,  EMG(1,M),  E.MG(2,M), 
EMG(3,M),  ENG(4,M)  Tl,  T2,  T3,  T4,  T5). 

The  arguments  Tl  through  T5  are  none.ssential  for  the  input  and  initial¬ 
ization  phase  and  are  set  =  0  before  the  CALL  to  POROUSl. 
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smamm  pags  buk-mov  num 


Table  B.l 


GLOSSARY  OF  MAIN  VARIABLES  IN  POROUSl* 


■I 

■| 


BKC  (K^) 

- 

Bulk  modulus  at  consolidation 

BKI  (Kj) 

- 

Initial  bulk  modulus  of  porous  material 

BULK  (K) 

- 

Linearly  varying  bulk  modulus  of  porous  material 

D  (p) 

- 

Current  density 

DC  (p^) 

- 

Consolidation  density 

DI  (Pj) 

- 

Initial  density  of  porous  material 

DIG  (p  ) 
so 

- 

Density  of  nonporous  material  at  ambient 
pressure 

BOLD 

- 

Density  at  start  of  current  solution  step 

LINT,  HOLD 

Current  value  of  internal  energy  and  value 
at  start  of  current  solution  step,  respectively; 
EOLD  is  not  being  used 

EMGl,  EMG2, 

EMG3,  EMG4 

Coefficients  of  Mie-Griineisen  equation  for 
nonporous  material  (EMGl  =  A,  EMG2  =  B  in 
Equation  6) 

ICON 

Array  of  consolidation  indicators  for  all 
the  cells.  ICON  is  0  before  consolidation 
and  1  after  consolidation 

JCELL 

- 

Cell  number 

M 

- 

Material  number 

>LATL 

- 

Array  that  identifies  material  m 

NFS 

- 

Number  of  parabolic  segments  defining  porous 

P-V  curve 

P 

- 

Current  pressure  computed  by  POROUSl 

PA 

- 

Pressure  in  the  porous  material  during  initial 
loading 

Items  in  parenthese  represent  the  corresponding  quantities  using 
the  notation  of  Section  3. 
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Table  B.l  (concluded) 


PB 

PMM  (p) 
PORA,  PORB, 

PI,  P2  (P^, 
SPV 

VI.  V2  (V., 


GLOSSARY  OF  MAIN  VARIABLES  IN  POROUSl 


Pressure  in  the  porous  material  during 
unloading  or  reloading 

-  Measure  of  volumeric  strain  (D  -  DIG) /DIG 

PORG  -  Goefficients  (constant,  linear,  and  quadratic 

terms,  respectively)  defining  the  parabolic 
segments  of  the  porous  P-V  curve 
PORA  =  P. 

PORB  -  -  Pi)/(Vj„  -  V.) 

PORC  =  -43P/(V.,,  -  V.)^ 

1+1  1 

P^^j^)  -  Pressures  at  the  endpoints  of  a  parabolic 

segment 

Array  of  specific  volumes  at  the  endpoints 
of  the  parabolic  segments 

-  Specific  volumes  at  the  endpoints  of  a 
parabolic  segment 
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Table  B.2 


FLOWCHART  OF  MAIN  EQUATIONS  IN  POROUSl 


(I)  time  =  0:  B  AD  in  data  and  compute  the  coefficients  of  the 
parabolic  segments  defining  the  porous  material  P-V  curve, 
i.e.,  for  parabolic  segment  i,  a,  b,  c  in  Equation  (3)  are 
given  by 


a  =  P. ,  b  = 
1 


P-^1  -  P- 
1+1  1 

V...  -  V.  ’ 
1+1  1 


c  = 


-4DP 


‘''i+i  -  “i' 


(II)  time  At: 

(1)  If  materia]  of  cell  j  i.as  consolidated  (ICON(j)  >  0  or 
D  >  DC),  go  to  (7). 

(2)  Determine  endpoint  specific  volumes  bounding  current 
specific  volume. 

(3)  Compute  pressure  for  initial  loading 

P  =  a  +  b  (V  -  Vj^)  +  c  (V  -  V.)(V  -  Vj^j)  (from  Equation  3) 

(4)  Compute  current  value  for  bulk  modulus  of  porous  material 

K  -  K, 
c  X 

K  =  K-  +  -  (p  -  p,.)  (from  Equation  5) 

I  Pj,  ~  ^I  I  ^ 

(5)  Compute  pressure  for  unloading  or  reloading  of  porous 
material 

pJi  _  pH  1  _  j,  ^  V/V  (from  Equation  4) 

avg 

(6)  Actual  pressure  in  porous  material  is  given  by  the 
minimum  of  (3)  and  (5).  Go  to  (8). 

(7)  Compute  pressure  for  consolidated  material, 

2 

P  =  Ap  +  (from  Equation  6) 

where  u  =  j/p  -1 
so 


(8)  RETURN. 


The  POROUSl  pressure  computations  occur  during  the  solution  portion 
of  the  code.  For  example,  SUBROUTINE  SWEEP  is  the  routine  in  COPS 
that  solves  the  governing  equations  (conservation  of  mass,  momentum, 
and  energy)  to  determine  the  particle  velocities,  displacement,  and 
strains.  It  transfers  control  to  the  various  material  models  to 
obtain  stresses.  Before  the  CALL  to  POROUSl,  NCALL  is  set  =1,  thus 
transferring  control  to  the  pressure  computation  portion  of  POROUSl. 

•A  listing  of  POROUSl  is  shoum  in  Figure  B.l. 
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FIGURE  B.1  LISTING  OF  SUBROUTINE  P0R0US1 


FIGURE  B.1  LISTING  OF  SUBROUTINE  POROUS1  (Continuad) 
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Appendix  C 
DEVIATOR  STRESS  MODEL 

The  basic  equations  for  the  stress  deviator  mode]  are  presented  in 
section  3.  Here  we  derive  the  relations  needed  to  use  these  equations 
in  the  numerical  calculations.  Listings  of  the  stress  deviator  sub- 
iwUtine,  VPSH,  and  the  subroutine  that  describes  the  strain-hardening 
yield  surface,  J2I2P,  are  presented  at  the  end  of  this  ap-  dix  along 
with  the  procedure  for  inserting  them  into  the  wave  propagation  code 
(COPS). 


ELASTIC  STRESS-STRAIN  RELATIONS 

In  the  elastic  region,  an  increine..tal  form  of  Hooke's  law  is  used 
to  determine  tiie  stresses 

•  da."  =  2Gde.:®  (C.l' 

ij  ij 


where 

cio.: 

ij 

=  da. . 
ij 

mm 

3 

6.  . 
ij 

de 

mm 

Qc  ,  . 

“  QC  .  . 
ij 

3 

0  ,  , 

IJ 

Here  G  is  the  shear  modulus,  assumed  to  be  a  function  of  pressure. 

PLASTIC  STRESS-STRAIN  RELATIONS 

A  yield  surface  of  the  type 

f(/jJ,  =0  (C.2) 
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is  assumed  to  separate  elastic  and  plastic  deformations.  In  our  model 
we  assume  plastic  incompressibility  and  isotropic  strain  hardening.  The 
yield  surface  consistent  with  these  assumptions  and  the  LLL  data  is 
given  by 


f  =  =  0 


(C.3) 


where  =  Yield  stress  expressed  in  terms  of  / 

M  =  nonlinear  strain-hardening  coefficient 


/Jo  =  /-r  O.  .  .  a . . 

2  2  ij  ’  j 

/If  -At.:’’  . 

2  2  ij  ij 

The  function,  f,  was  obtained  from  the  uniaxial  stress  data  obtained  by 
Heard  et  al.  The  strain-hardening  plot  (/^J^  ~  versus  / obtained 
from  their  data  is  shown  in  Figure  C.l.  This  curve  was  fitted  by  a 
series  of  parabolic  segments  similar  to  those  used  in  the  porous 
pressure-volume  relation  (Section  3).  The  parameters  describing  these 
parabolic  segments  are  given  in  Table  C.l.  The  strain-hardening 
coefficient,  M,  for  a  particular  strain  increment  can  then  be  obtained 
from  the  secant  modulus  appropriately  centered  over  the  increment. 

For  a  proportional  loading  experiment  (like  uniaxial  stress),  the 
plasticity  model  used  here  gives 


(C.4) 


where 


'"ai^  =  /j  AeT?  .  AeT? 

2  2  ij  ij 


Equation  (C.A),  derived  in  the  next  subsection  of  this  Appendix,  is 
required  for  the  strain-hardening  model  used  here. 


H.  C.  Heard  et  al.  "Stress  Strain  Behavior  of  Polycyrstalline  NaCl 
to  3.2  GPa,"  Lawrence  Livermore  Laboraotry,  Livermore,  C-^lif., 
Report  No.  HCRL-51743  (1975). 


Table  C.l 


PARAMETERS  DEFINING  THE  SIX-SEGMENT  PARABOLIC 
APPROXIMATION  TO  THE  STRAIN-HARDENING  CURVE 


yJZ  -  Y  (kbar)  A[/j^  -  Y  ]  (kbar) 

2E  q _  _ — - - - 


0.05197 

0.10393 

0.195 


0.218 


0.007 

0.01366 


0.004 


0.000965 

0.0042395 

0.021223 


0.037882 


0.068 


Y,^1  is  the  maximum  deviation  from  a  linear  "  Y^  versus 


'ly  variation. 


At  the  start  of  a  computational  cycle,  the  deviatoric  stresses 
are  formed  by  assuming  an  elastic  increment 


a.%  =  a..'  +  2GAe.: 

ijl  ijo  13 


(C.5) 


where 


a.'.  =  stress  deviator  at  the  end  of  previous 
i  10  * 

computational  cycle 

G  =  current  value  of  shear  modulus 


Equation  (C.5)  is  used  to  compute  /j'^.  If  / is  not  greater  than 
the  current  yield  strength  / (or  Y) ,  no  adjustment  is  made  to  the 
elastic  stress  values  determined  from  (C.5),  and  the  calculation  proceeds 
I£/J^  is  outside  the  current  yield  surface,  stresses  are  determined 
from  the  relation 


-  n  "  1  /i  \At 

j  =  "iji  I' -73^/^ 


(C.6) 


where 


To  =  relaxation  time 
K 

At  =  time  increment  <  T 

K 


For  no  strain-hardening,  all  quantities  on  the  right  side  of  (C.6)  are 
known  and  a. .  can  be  evaluated.  Plastic  strain  increments  are  obtained 
from 


a.  -  a. .  a7._ 

p  =  ^3 1  1.1  ial 

j  2G  2G 


he..  = 

ij 


where 


a.  .  =  o.  -rrp 

ij  ijl  V J. 


:4- 


In  the  pre.>ence  of  strain-hardening,  the  value  of  / used  in 
(C.6)  must  be  consistent  with  the  plastic  strain  generated  during  the 
increment.  Strain  liardening  during  a  cycle  is  expressed  as 


=  °/j^^  +  M  /ai;’’ 


(C.8) 


where  V is  the  yield  strength  at  the  start  of  the  current  increment 
and  M  is  the  current  work  hardening  modulus-  Using  (C.7),  we  can 
x%'rite 


2G 


(C.9) 


Equation  (C.6)  can  be  expressed  as 


/jr  = 


/j'  (1  -  +  /.I 


1' 


2E 


At 

T, 


R 


(C.IO) 


Equations  (C.8)  through  (C.IO)  have  three  unknowns, 

Hence  a  simulrnneous  solution  can  be  obtained  to  give 


/  t'  o 

'  ^21  -  '  "2E 
.  2GT„ 


At 


(C.ll) 


Substituting  (C.ll)  in  (C.8),  we  get  This  value  when  substituted 

in  (C.6)  gives  o.'..  Tiie  plastic  strain  increment  AcT^  can  then  be 
ij  ij 

determined  from  Equation  (C.7),  and  the  calculation  proceeds. 

The  simultaneous  solution  presented  here  is  possible  because  of 
Equations  (C.4)  and  (C.8)  and  the  determination  of  M  from  Figure  C.l. 
The  last  two  items  are  discussed  here  and  the  derivation  of  Equation 
(C.A)  is  given  in  the  next  subsection  of  this  appendix. 

_  /~Ar 

M  is  a  nonlinear  function  of  *'1^  . 
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. . .  ^ 


! 


Equation  (C.8)  represents  a  linear  work-hardening  relation  over 
the  strain  increment.  This  linear  relationship  is  justified  because 
we  are  using  the  current  value  of  M,  because  strain  hardening  is 
monotonically  increasing,  and  because  we  are  using  an  incremental 
approach. 

The  determination  of  M  from  Figure  C.l  relies  on  Equation  (C.4) 
and  the  approximation  AeT*?  ~  ^^ij*  it  is  possible  to  improve  on 

this  approximation  by  iterating  on  M  using  the  value  obtained  from 
Equation  (C.ll),  the  difference  in  the  results  is  insignificant  because 
of  the  relative  magnitudes  of  M  and  G  in  Equation  (C.ll). 


DERIVATION  OF  EQUATION  (C.4) 


=  A 


in  _  , 


1/2 


(C.12) 


The  material  model  considered  here  is  equivalent  to  the  Prandtl-Reuss 
relations,  which  give 


de'^  =  aC..dX 

iJ  ij 


(C.13) 


where  dX  is  a  nonnegative  scalar.  Substituting  (C.13)  on  the  right 

A 

side  of  (C.12),  we  have 


For  proportional  loading,  we  can  write 


(C.14) 


a'.  =  6..0  (C.15) 


*  n  * 

Note,  de*^  =  0. 
mm 
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. . . 


where  o  a  proportionality  constant.  By 

taking  the  dot  product  of  (C.15)  with  itself,  we  can  write 


6..  -  6..  =  1 

2  ij  ij 


(C.16) 


Substituting  (C.15)  in  (C.14),  we  can  take  out  3..  from  the  resulting 

^  J 

integral.  Using  (C.16)  we  can  write 


Av/l^  =  Ai/o  dX] 


(C.17) 


If  d(odX)  is  defined  over  the  limits  of  integration,  the  right  side  of 
(C.17)  is  simply  0  dX.  Using  the  definition  of  o,  we  have 


=  adX  = 


k  o:.dx  .  o:.  dx 

2  IJ  IJ 


(C.18) 


Using  (C.13),  we  can  write 


A/l^ 


lA  Ae:.'’  .  Ae:.P  = 

»  2  IJ  n  2 


(C.19) 


The  equivalence  shown  in  (C.19)  is  possible  because  of  the  proportional 
loading  assumed  in  (C.15).  While  the  compression  and  shear  experiments 
are  nonproportioual  loading,  the  work-hardening  function  is  constructed 
from  the  static  uniaxial  stress  experiment,  which  involves  proportional 
loading. 


VPSH  AND  J2I2P  INSERTION  PROCEDURE 

VPSH  is  a  two-dimensional  visco  plastic,  strain-hardening  model. 
J2I2P  defnines  the  isotropic  strain-hardening  curve  used  by  VPSH.  Both 
subroutines  consist  of  two  section,  a  data  input  and  parameter 
initialization  section  and  a  coraputat'on  section. 


m 

r  -trJ 


hJi 


The  argument  lists  of  VPSH  and  J2I2P  are  as  follows; 

SUBROUTINE  VPSh  (NCALL,  IN,  J,  M,  MATE,  YO,  Yl,  Y2,  GO,  Gl,  G2,  DT,  SXXT, 
SYYT,SZZT,  SXYT,  EXP,  EYP,  EZP,  EXYP,  EX,  EY,  EZ,  EXY, 
INDY,  DRUG) 

and 

SUBROUTINE  J2I2P  (NCALL,  R,  J,  M,  WHMO,  I2P,  J2E,  DRUG). 

Most  of  these  parameters  are  defined  in  Tables  C.2  and  C.3.  Table  C.A 
outlines  how  the  equations  presented  above  and  in  Section  3  are  used  in 
the  numerical  calculations. 

In  both  subroutines,  NCALL  directs  control  to  either  the  data 
input  and  parameter  initialization  section  or  the  computation  section. 

For  NCALL  =  0,  VPSH  reads  in  the  required  data  from  input  film  IN 
(usually  =  5)  and  CALLS  J212P,  which  in  turn  reads  in  the  data  that  it 
requires  from  input  file  K(=  IN)  and  initializes  the  necessary  parameters, 
i.e.,  the  coefficients  defining  the  parabolic  segments  that  approximate 
the  isotropic  strain-hardening  curve,  and  the  initial  work-hardening 
modulus.  For  NCALL  =  1,  VPSH  computes  the  relaxed  deviatoric  stresses 
and  plastic  strain  increments.  In  the  process  of  doing  so,  VPSH  CALLS 
J2I2P  with  values  of  to  obtain  corresponding  values  of  ^ 

VPSH  and  J212P  data  input  and  parameter  initialization  are  performed 
as  part  of  the  problem-generating  portion  of  the  code  (e.g.,  SUBROUTINE 
SETUP  in  the  COPS  code).  Before  the  CALL  to  VPSH,  NCALL  and  IN  are 
set  equal  to  0  and  5,  respectively,  indicating  that  control  is  to  be 
transferred  to  the  data  input  and  parameter  initialization  sections  of 
VPSH  and  J2I2P.  The  CALL  to  VPSH  in  SETUP  is 

CALL  VPSH  (NCALL,  IN,  J,  M,  MATL,  YO(M) ,  Y1(M),  Y2(M),  GO(M) ,  G1(M), 

G2(MI,  Tl,  T2,  T3,  T4,  T5,  T6,  T7,  T8,  T9,  TIO,  Til,  T12, 

T13,  TIA,  IDBUG). 

The  arguments  Tl  through  TIA  are  nonessential  for  the  input  and 
initialization  phase  and  are  set  equal  to  0  before  the  CALL  to  VPSH. 


The  VPSH  computations  of  deviatoric  stresses  and  deviatoric  plastic 
strain  increments  occur  during  the  solution  portion  of  tne  code  (e.g., 
SUBROUTINE  SWEEP  in  the  COPS  code).  VPSH  is  called  in  the  portion  of 
SWEEP  that  treats  the  c'eviator  stress  material  models.  Before  the  CALL 
to  V'PSH,  NCALL  is  set  equal  to  1,  thus  transferring  control  to  the 
computation  sections  of  VPSH  and  J2I2P. 

Figure  C.2  shows  a  listing  for  subroutine  VPSH,  and  Figure  C.3 
shows  a  listing  for  J2I2P. 


. . . 


i  1 


DEV  (-y51) 


Table  C.2 


GLOSSARY  OF  MAIN  VARIABLES  IN  VPSH 


-  Volumetric  strain  increment 


DEX,  DEY,  DEZ,  DEXY 


(dc'^,  de'^) 

X  y  2  xy 

nT9D  /'yAT'F^ 


DI2P  (.'AI'*  ) 


DT  (At) 


EX,  EY,  EZ,  EXT 


(de  ,  de  ,  dc  ,  de  ) 

X  y’  z  xy 

EXP,  EYP,  EZP,  EXYP 

(l'P,  f'P,  C' 

X  ’  y  ’  xy 

EXPP,  EYPP,  EXPP,  EXT’PP 


GO,  Gi 


I2P  ot;P) 


-  Plastic  deviatoric  strain  increments 


Second  invariant  of  the  plastic 
deviatoric  strain  increments 


I2PO  (°/l'”)  • 


-  Current  time  step 


-  Total  strain  increments 


-  Plastic  deviatoric  strains 


-  Temporary  plastic  deviatoric  strains 
obtained  by  assuming  that  the 
dev'iatoric  strain  increments  are  all 
plastic 

-  Constant  and  linear  shear  moJulu- 
coef ficients 

-  Two  times  the  shear  modulus 

-  Second  invariant  of  the  plastic 
deviate. ic  strains 

-  Second  invariant  of  the  temporary 
plastic  deviatoric  strains 

-  Second  invariant  of  the  plastic 
deviatoric  strains  at  the  start  of 
the  current  cycle. 


s 

Items  in  parentheses  represent  the  corresponding  qua.itities  using  the 
notation  of  Section  3  and  r’^is  Appendix. 


i  I 


I 


i  * 

I  I  U 


fm 


Table  C-2  (continued) 


GLOSSARY  OF  MMX  VARIABLES  IX  VFSH 


J2E  (/J'  ) 


J2E0  (Vj;^) 


J2T  (vOT) 


RELAXF 


-  Second  invariant  of  the  equilibrium 
deviatoric  stresses  (computed  from 
Kq-  C.8) 

-  Second  invariant  of  the  equil ibrium 
deviatoric  stresses  (corresponding 
to  I2P) 

-  Second  invariant  of  the  equilibrium 
deviatoric  stresses  corresponding  to 
the  temporary  plastic  deviatoric 
strains 

-  Second  invariant  of  the  equilibrium 
deviatoric  stresses  at  the  start 

of  the  current  vcle  (corresponds  to 
I2P0) 

-  Second  invariant  of  the  deviatoric 
stresses  computed  assuming  all 
elastic  strain  increments 

-  Material  number 

-  Array  that  identifies  material  m. 
Relaxation  function 


L-  (l-  73^ 


SXX,  SYY,  szz,  5XT 


-  Relaxed  deviatoric  stresses 


k 

Items  in  parentheses  represent  the  corresponding  quantities  using  the 
notation  of  Section  3  and  this  Appendix, 
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MM 


Table  C-2  (concluded) 


i 


GLOSSARY  OF  MAIN  VARIABLES  IN  WSH 


SXXT,  SYYT,  SZZT,  SXYT 

°>'r  %-i> 


TAUO  (T^) 
WHM  (H) 
WKMO 


-  Deviatoric  stresses  computed 
assuming  all  elastic  strain 
increments 

-  Set  equal  to  SXX,  SYY,  SZZ  and  SXY 
before  RETURN  from  VPSH 

-  Relaxation  time  constant 

-  Work-hardening  modulus 

-  Work-hardening  modulus  at  start  of 
current  cycle 


•  I 


* 

Items  in  parentheses  represent  the  corresponding  quantities  using  the 
notation  of  Section  3  and  this  Appendix. 
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Table  C.3 


GLOSSARY  OF  MAiX  VARIABI£S  IN  J2l2r 


Cl,  C2.  C3 


-  Coefficients  of  constant,  linear, 
and  quadratic  tenas,  respectively, 
of  parabolic  segaents  used  to 
approxinate  the  strain- hardening 


C!  . 


DV  -  v^l.) 


-  Maxiaun  deviation  froa  a  linear 

yi  ■!  '■ 

(/Jl_  -  Y  )  versus  /I-'  variation 
,  2E  o^  .  2 

for  sccaent  i 


I2P  (.'i;") 


-  Second  invariant  of  the  deviatoric 
nlastic  strain 


J2E 


-  Second  invariant  of  the  equilibrixrcs 
deviatoric  stresses 

-  Corresponds  to  12P 


-  Material  nuaber 


YO  (Y  ) 
o 


-  Nunber  of  parabolic  segaents 

-  Initial  work-hardening  ssodulus 


-  Array  containing  the  endpoint 

values  of  -  Y' 

2t  o 

-  Array  containing  the  endpoint 
value  of 

-  3.t  initial  yield 


Iteas  in  parentheses  represent  the  corresponding  quantities  using  the 
notation  of  Section  3  and  this  appendix. 


Table  C.A 


FLOWCHART  OF  MAIN  EQUATIONS  IN  VPSH  AND  J2I2P 


(I)  time  =  0:  READ  in  data  in  VPSH,  CALL  J2I2P  to  READ  in  data  and 
compute  the  coefficients  of  the  parabolic  segments  defining  the 
strain-hardening  curve,  i.e.,  for  segment  i,  a,  b,  c  in  Eq.  (3) 
are  given  by 


=  (/J'  ),,  b  = 


^'^‘^2E^i+l  '  ^*^'^2E^i 


2E  1 


^  ^2  ■'i+1  ^*^^2  -'i 


c  = 


<■'12  >i«  - 


(II)  time  >  At: 


(1)  Compute  the  trial  deviatoric  stresses 


a.'  =  a:.  +  2G  A  eC. 

ijl  ijo  IJ 

(C.5 

(2) 

Fora/J^. 

(3) 

It  /Tj;  i  VT;;  -  4  e-/  .  0  * 

go  to  (15) 

(M 

Form  Vl^  =  -l/i  e:.P  e-:.P 

L  »  2  ijo  ijo 

(C.3 

(5) 

CALL  J2I2P  to  obtain  i' 

IE 

(6) 

Compute  the  plastic  deviatoric 

strains  assuming  the 

deviatoric  strain  increments  are  all  plastic 

=  c'.  y  +  AeT. 

IJO  IJ 

(7) 

Form  /l^^  =  y^4 

2  »  2  ij  ij 

*This  is  done  in  the  COPS  SUBROUTINE  SWEEP. 


Table  C.4  (continued) 


FLOWCHART  OF  MAIN  EQUATIONS  IN  VPSH  AND  J2I2P 


(8)  CAI.L  J2I2P  to  obtain  /j: 


(9)  Determine  the  work-hardening  modulus 

M  =  vj^)/(Vi^P  -  vTJP) 


IF(  /l^P  -  Vl^)  <  0,  M  =  “M 


(10)  Compute 

_  _  _  2GT 

/M'P  -  (/J'j  -  Vj'^)/<M  +  -^) 

(11)  Compute  /j^  =  Vj^  +  M 


(12)  Form  the  relaxation  function  and  compute  the  relaxed 
deviatoric  stresses 

r  /J^l 

^  1  2E  At-. 

0 . .  =  0 .  _  {1-  1-  -7-p- 

xj  ijl  I  J  \ 


(13)  Compute  the  plastic  devia*"'rlc  strain  increments 
and  update  the  plastic  deviatoric  strains. 


(C.ll) 


(C.8) 


(C.6) 


Ae:P  =  (a:.^  -  a:.)/2G 
ij  ijl 


eT?  =  eT.P  +  AeT^ 

ij  ijo  IJ 


(from  Equation  (15)  and  (C.7) 


This  is  done  in  the  COPS  SUBROUTINE  SWEEP. 


■1 


H 

i'l 


I 


■4i 


Table  C.4  (concluded) 

FLOWCHART  OF  MAIN  EQUATIONS  IN  VPSH  AND  J2I2P 


(14)  Update  the  work-hardening  modulus  for  the  next  cycle 


(14. a)  Form /T?  =  eC?  eC? 

2  »  2  13  ij 


(14. b)  CALL  J2I2P  to  obtain  /j^T 

2E 


(14.  c)  °M  =  (/J^g  -  V  ^ 


If  -  Vl^)  <  0,  °M,  =  M  from  (9) 


(15)  RETURN 
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fO- 


subroutine  VPSH  COC  6700  FTN  V3.0*3S5F  OPTcl  7’>/06/a5 

SZZTbSZZ 
SX7  T=SXY 


FIGURE  C.2  LISTING  OF  SUBROUTINE  VPSH  (Concluded) 


FIGURE  C.3  LISTING  OF  SUBROUTINE  J2I2P  (Concluded) 
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